ENU 6052 - RADIATION TRANSPORT BASICS & APPLICATIONS 


INTRODUCTION 


Statistical Nature of Matter 


Typical macroscopic descriptions of physical phenomena are 


N(t) = ett 


i.e., the population of radioactive nuclei remaining at time t, and 


I(x) = te =* 


i.e., the intensity of a monovelocity particle beam reamining uncollided to 

a penetration x in an interactive medium. There are several conceptual diffi- 
culties inherent in these results. To resolve some of them, and provide an 
introduction to the material of these lectures, consider the radioactivity 
problem in more detail. 

"As good as new" postulate: 

The probability that a radionuclide decays in the time period t to t + At 
is \dt where A is a constant which depends on the nuclear specie. Note that 
At must be "small," i.e., AAt << 1 and that the probability is t~independent. 
The quantity of interest, in studying the detailed meaning of the N(t) ex- 
pression, is e (At,t) = probability. that m nuclei decay in t. to * + At pro~ 


viding that there are % nuclei present at te. i NeEe that, for ‘smal: At, _ che 


"as good as new" postulate for each nucleus ‘implies EERE 2 Maes 
Binomial distribution: 
be) = probability. that m entities (out of Doca fy 

—& is the probability that each experience the event. 7. ae 


L Xt ~ 
be) = argemT Sma-ey™ 
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and has the properties 


& 
dy bd (0) = Sn0 
& “ 
2. biG) = Of 
. fly 455. 
where 5,5 {i ij 


£8 
3. Eid (6) = 1 


tr ok 
4, % pmb (6) = 28 


5. 2p (e) = RE + 2(g21) Ee 
Ei gmb ice) = £(221) E 


Examples are A(— = 0.5) Océ = 0.8) 
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#0 =—50 


reds 


Poisson distribution and normal distribution: 
For € << 1 and 2 >> 1, 


m -2& 
bce) > pace) = Ge 


mi 
which is the Poisson probability distribution. An alternative approximation 
is the normal distribution which is defined for the continuous variable pu, 


i.e., 


discrete m + continuous yp, 


exp E(u 28) /2.07 


g 
p.(é) = 
u To 


where o* = 2&(1-£). pice) is a continuous probability distribution in the 
sense that 
pC edu = probability that y» in the range p to p+ dp entities out of 2 


experience the event. 


Example comparison of b(), pate) and pA (é)- f= 


Radioactivity relations: 


at(at) = bicade) for AA << 1. 


Consider N, radionuclei present at t = 0 and define 


Pi (t) 
N 


probability that j nuclei have not decayed by time t. 


lo 
N(t) = F20 JP, (0) = "expected" (i.e., average) number of nuclei which 


have not decayed by t. Question - What temporal relation does N(t) satisfy? 
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Solution - For small At, 


j At) 
LqP.(e + At) = LD 4p. (ryqitm ¢ 
43 act t) Fee jm! ) 


& 2 
= Roney 2 
= aed ozs mo %,(A®) Be, (t) map BQ, (At) 
z 2 
z = 
But, m0 Qa at) 1 and 
g Q n L be be 
mo MQ, (At) = 29 mb CAdt) = grat 
whence, 


N(t + At) - N(t) = -An(t)At 


N(t_+ At) = N(t) 


At = ~AN(t) 


and from the limit as At+0, 


aN (t) 
oe Coa al 


with solution N(t) = wie 


In the familiar radioactivity relation, N(t) actually means the statis- 
tical average result obtained with many measurements commencing with identical 
initial conditions. 


Particle transport problem: 


7 ‘ 
| 

4 

pomenenn ey , ‘ 
' 

= 
4 

ney ¢ = 


oN 
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Consider a beam of monovelocity particles perpendicularly incident on a 
half-space of targets of uniform density. For each particle, the probability 
that the first interaction with a target occurs between x and x + Ax is x- 
independent and, if Ax is small enough, is given by ZAx. If (bse = proba- 
bility that m particles have first interaction in x to x + Ax with & uninter- 
acted particles at x, then (Axe) = ba (Ede) And if I(x) = FIP; GOs then 


wee = ~DI(x) implies I(x) = Le, 


Coordinate Transformations 

Vector x means (x. 2Xpreeer%)> e.g., location on a line x means x, 
Cartesian location on a plane x means (x,y), Cartesian location in configura- 
tion x means (x,y,z). A vector x can be expressed in alternate coordinates, 


i.e., x can be represented by 


Cy 2Xq 90+ +9) or (a ,Uyseeest)s 


e.g., location on a plane. 


vr r=r (x,y) 
(x,y) 

we = 8 (x,y) 
e (1,0) 


and if angle limits are set to 
yield one-to-one transformation 


x= x (r,6) y = y (1,6) 


Differential element d"x means dx, dx dx,...dx,, @.g-, line aby means 

dx, Cartesian plane dx means dxdy, polar plane d*x means drd8, spherical 
polar configuration d*x means drdOdd. The function f(x) means E(x Xp, -+ ek } 
= n 


and if the function is itself a vector, £(x) means [£, (xy ---x,)5 oes £ (Rye KT 


Density Functions: 

f£(x) will often represent the density of some physical property in the 
space of x. 
Examples are 

1. Mass density along a line, i.e., f(x) such that the mass between x=a 
and x=b (a<b) is given by P £(x)dx. 

2. Particle density = configuration, i.e., f(x,y,z) such that the number 
of particles in volume V is given by f £(x,y,z)dxdydz. 

3. Particle kinetic energy spectrum, i.e., f(£) such that the number of 


particles with kinetic energy between E, and Ey (Ey < E,) is given by 


1 


f"2 eceyan. 
Ey 


Transformation of coordinates: 

One-dimensional case - 

Let u = u(x) be a one-to-one transformation (x*u) with inverse x=x(u) 
describing transformation (u»x), and f(x) be a density function such that 
Fe ;° £(x)dx has some physical significance. Problem - To express F as an 
ietecres over u rather than x and thereby deduce the form of the density 
function in u-space. 

Solution ~ x = x(u) implies dx = au du 
and 


(b) dx 
Fe f°? gtx(uy] & au 
ues du 


However, it is customary to choose integral limits such that 


(lower limit) < (upper limit). 


Moreover, if u(x) is one-to-one, then it is monotonic. Whence, 


f= fM*" elx(u)] el du 
min u 
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and conclude that 
g(u) = £[x(u)] |S 
is the density function in u-space such that F = f g(u) du, e.g., E = lonv? 


(v 2 0) and £(v) represent particle speed spectrum. 


Vv 

F= | - f(v) dv = number of particles with speed in range vy to v 
V1 
The corresponding particle kinetic energy spectrum is given by 


ae) - ¢ (PS (e- 


2° 


and 


E 
2 
Fe | g(E)dE, Ey = isnv? 


An explanation of the varying form of spectrum (density function) is aided 
by noting that equal intervals of v are transformed into unequal intervals of 


E as illustrated. 


LH 
ZO 
/6 


nw Ee 
a 
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Thus, for example, it is expected that a constant density (spectrum) in 
speed space implies an energy spectrum which decreases with increasing E, 
Two-dimensional case - 
Let u = u(x,y), v = v(x,y) be a one-to-one transformation (x,y) + (u,v) 


with inverse x = x(u,v), y = y(u,v), i.e., (u,v) + (x,y). 


vy 
sf] 


[f £Gx,y) dxdy 
R 


hy 
u 


H £[x(u,v), y(u,v)] [2X fauav 


where the Jacobian of the transformation is defined by 


ox Bx 
KY: gu ov 
oe 2 
du ov 


The density function in (u,v) - space is 


g(u,v) = f[x(u,v), y(u,v) } [52% | 


Q.Bes 


Thus, g(r,9) = f(rcos6, rsing) r 


r= [f g(x,8) drde, or equivalently, 
R 
Fe ff f(reos6, rsin6) r drdé 
Se ea SORRY. eT ee 


R * * 
*Note that the indicated grouping of factors is not in the spirit of ex- 
pressing density functions and differential elements in (r,6) ~ space. Rather, 


the grouping implies merely a change of coordinates for performing the integra~ 


tion ff f(x,y) dxdy, and 
R 


x dr dO is an area differential of the 


type dxdy, i.e., in (x,y) - space, not 


(r,8) - space. 


n-dimensional case - 


Let x = x(u) have inverse u = u(X). 
F= f £(x)d"x = f g(u)d@u 
R R 


where 


g(u) = £[x(u)] | J} 


Ox, 8x, | | OR, 
2 


ou, du du, 
x ax 
I@ = 2 
ae Lon 
iL, 
ox, 3x, 
Ta, ei Beoe da 
n 


Particle Density Distributions 

Configuration-space density: 

Consider N particles labelled g = 1, 2, ..., N. The’configuration of 
these particles in 3-space is completely specified by the set of 3N numbers 
{r(o), a=1, 2, ..., N}. 

There are at least two reasons why such a dynamic state description is 
unrealistic - 

1. N of interest is often larger than 107°, 

2. Particle interactions are thought to be actually statistical in 
nature (quantum mechanics) and thus a precise specification is of 
questionable relevance. 

Thus, there is motivation for developing the ideas of describing particles 

by statistical distributions and their moments (e.g., averages). For example, 


the one-dimensional particle configuration 


6 —e- o> Xx 


@ 2 @~ 
x(t) x) x (3) x(4 x08) x6) x7) 


is presumably unpredictable (i.e., with given initial conditions, identical 
experiments yield different results at time of measurement). Instead of 
specifying the set {x(a)} a particle density function is developed and defined 


in a probabilistic fashion. 


For example, suppose the following is obtained in 5 measurements. 


4 | nce Gh omit ae 
| 1 
; 4 
. Z bitten ate 
NY | i 
; | 
~ ; i : j | i 
& 3 -@-8—-2-01 9-6-2 ---@—~9 Le 
oO me ee 9 _—_ 
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Nt i ' i 
ls L : { : 1 | i i 
g pe tent fener te penne 
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N i i | : 
| 
s le-@—@--ei—@-- 6 pe Peete eee Be } e oo ee “ 
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o Ax ZAK FAY ¥Ax SAX é4x 7AX 


Avg. KesulA 3,2 22 ps4 he ho 0,8 26 


With respect to this set of measurements, define 
Ry (x,A4x) = probability of finding j particles in Ax at x. 


If there are many more than 5 measurements performed, then often 
lim P.(x,Ax) = S,, 
dxro 4 i 


and the limit is approached smoothly such that 


N(x) = lim 1 24P (x, Ax) 
Ax +0 Ax jj 


has a clear meaning, 


b « 
f N(x)dx = N(a<x<b) = "expected" (or, average) number of particles 
a 
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found in the interval a © x § b over many identical performances of the experi- 
ment. Generalizing to n-dimensional configuration space (n=l, a line; n=2, 


a surface; n=3, a volume). N(r) is so defined such that 
Syqeyae = N(R) 
R 


i.e., the expected number of particles in region R. 


Transformation of configuration density: 


Let r' = r'(x) be a one-to one transformation with inverse r = r(r'). 
Then, 
wen) = [w@ ate =f nt@att 
R R 
where 


N(x!) = NG) |I@| 


Note that N' represents the particle density function in x' - coordinates 
and is usually a different function of its variables than N is of its vari- 
ables. 

For example, consider the 3-space transformation from Cartesian coor- 


dinates (x,y,z) to spherical polar coordinates (r,6,). 


= 
x = rsinécos¢d 
2) x f y = rsinésing 
z = rcos@ 
~ y oSeSn 
Pht 2 0S¢S2m 


X,¥.Z.) _ 2 8 
Ie rsin 
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N'(r,6,$) = N(rsin€cos$, rsin§sind, rceos6)r"sing where N'(r,6,¢) is the 
particle density in spherical polar coordinates, and N({x,y,z) is the particle 
density in Cartesian coordinates. 

Moments of the configuration density: 

Zeroth-order moments — 

Restricting attention to the usual 3~space density N(r) with r = (fy o¥y0%q)> 
define the zeroth-order moments as 


ne fff N(x) 5To5Tq) drjdrjdr, 


N, (x) = ff N(x) ,r),.r5)drodr, 


and similar definitions for 


N, (fy) and N,(r3), and 
Ny (typ) = [ N(x, ,r9,7,)dr, 


and similar definitions for Ny (9503) and Ny (e3r4)- Note that if a region 
(limits) are unspecified on an integral, then the usual meaning (in these notes) 
is to cover the entire relevant space of the variable. For example, the 


Cartesian 3-space zeroth-order particle densities are 


+o +00 +00 
N =f JL ek N(x, y, 2) dxdydz 


and represents the total number of particles (everywhere). 
400 +400 
N (2) = ff NGx,y,z)dxdy 


and represents the particle density per illustrated relevant volume element 
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foo 
Ny(ysz) = f° N(x,y,2)dx 


and represents the particle density per illustrated relevant volume element. 


As another example, consider spherical polar 3-space and designate the particle 


density in spherical polar coordinates as N'(r,6,¢). 


n= f° [" 7" (x,0,0) dradae 
a t) i) 


and, again represents the total number of particles. 
T p20 
Nic) = ff N'(R,0,4) dodg 
0 0 


and represents the particle density per illustrated spherical shell volume ele- 


ment. 
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Note that if the Cartesian 3-space density is a constant, i.e., 
N(x,y,z) = C, then 
N'(r,0,6) = Crsind 

and 
N, @) = cr? f a sin@d0do 


4n rc 


i] 


Higher-order moments - 
Recall the one-dimensional example of particles distributed:along the 
x-axis. A particular measurement yields {x(a), a=1, ..., N} and the 
familiar definition of the average value of displacement x is 
N 
= 


1 
<x = = 
* N a= 


1 x(a) 


For the equivalent particle expected distribution, N(x)dx = number of parti- 


cles in dx about x and 


1 oo 2 
<x> = rf oy xN (x) dx 


Generalizing (in one-dimension), 
L co 
<a(x)> =H f eC) N(x)dx 
—_ 
where g(x) is any meaningful description of the particle location. 


Generalizing to n-dimensions, 


<g(z)> = Jg()N@ae 


and, if g is a vector quantity itself, 
<g(e)> = = faQN(ed"r 
Examples — 


1. Let g(r) = r in 3-space and use spherical polar coordinates. 


<> = f J" Jun’ (x, 0, $)drdedd 
a) 0 


or, in terms of the Cartesian density N(x, y,z) 


vo 4. 277 
<r> = é fff IN(x, y, z)x°sinSdrd8dd 
00 (0 


or, equivalently 


+o 400 too 2 


(x? +y + 27) N(x, y,z) dxdydz 


<x= & f° fT f?T ent (x, 0, 6)drasde 
0 0 0 


4n "isotropic" distribution is defined by N = N(x). For an isotropic 
distribution <r> = 0 but <r> = 0 only if all particles are at the origin of 


coordinates, 
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Configuration-velocity space density: 

For each particle (labeled @) r(@) = Zy(t) and v(@) = Vy Ct) are the 
particle trajectory relations and, thus, r(@) is a function of v(@). In 
order to describe the expected distribution of particles in both configuration 
(location) and velocity (energy and direction-of-travel), consider the 6- 
dimensional (r, v) ~- space where r = (> Xo> r5) describes location and v 
= (v4 5V52V3) describes velocity. Note that in this abstract space r and v 
are considered as independent variables. The particle density n(r,v,t) is 
defined such that [ ‘ n(xr,v,t) @raty = the expected particle population 
in configuration seid YY, with velocities in the range of volume vy in v-space. 
An alternative particle distribution description is in terms of zr, E and 8 where 
E is the particle kinetic energy and g unit vector in the direction-of-travel, 


Maw 
i.e., E= Semv? and 2 = 7, where m is the particle mass. The particle density 


n(x, E, &, t) is then defined such that 


E na 
f ae n(z,E, 0, t)d-rdka"2 
yk © 


is the expected particle population in configuration volume V with energy in 
the internal Ez to E, and with direction-of-travel in the solid angle W, at 


time t. 
Examples of 6-dimensional coordinate systems - 


1. (2,v) - space, Cartesian in both ¥ and v. 


%3 
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Atv 


4% 


AW 


2. (x,E,2) - space, Cartesian in r and spherical polar ink, 


IT 


£ 
| 
x 
AE 
dxdy az 


Transformation of (r,v) coordinates: 
Let r' = r'(r,v), v' = v'(r,v) be a one-to-one transformation with in- 


verse r = r(c',v'), v=v(r',v'). Then a'(r',v',t) = n[x(x’,v'), v(c',v"),t] 


EX. 
lI | A usual special case is the simpler transformation =r’ (r), 
> eon 


x 
vi =y'(y). Then n'(x'yv'5t)= afe(e'), viv'),t] [Gol] IG | 


Moments of n(r,v,t): 


N(z,t) = fn ny, t) dv = number density in configuration space. 
N(v,t) = fac, t)ar = averaged velocity spectrum. 


N(t) = Jface,v, t)d?ra2v = total particle population. 


g(x, t) = wee fe(eyy, t)n(x,v,t) dv are v-space moments (or averages) of 
is 


relevant physical quantities g. 


jo 
aaa 


<p (t) = ey ffeG@.y.t) rd“v are (r,v)-space moments (or averages) 


over total particle populations. 
Delta Functions 
Discrete index delta function: 
San has previously been defined in these lectures as 


§ = 1 if m=n 
mm 0 if mfn 


where m and n are integers and are thus a discrete set. An important feature 


of San is the collapse of summations to a single term, i.e., 


Note that if the set of terms {a} is arbitrary, then this collapsing property 
leads directly to the definition and can therefore be itself considered a 
definition. 

Dirac delta function: 

The Dirac delta function plays the same role with continuous variables 
that rai plays with a discrete index. In one-dimension, §(x-a) is defined 


by the relation: 


X, 


x 


7 £ 
f[ S(x-a)£(x)dx = | 
1 


(a) if X $a < % 


QO if a<x, or a>x 


1 


2 


which implies that §(x-a) can be considered as the limit of an "ordinary" 


function F(x,a) with the properties 


1. 


2. 


{° F(x,a)dx = 1 


oo 


In the limit, F(x,a) = 0 everywhere in x except near x 


it must be highly peaked. 


In n-dimensional space, §(x-a) is defined by 


[8(x-a)£(x)d"x = 
R 


f(a) if ainR 
QO if ais not inR 


For example, in Cartesian 3-space, the definition is satisfied if 


6(r-a) = 5(x-a) 6(y-a,) 6(z-a,) 


a where 


Dirac delta functions are density distributions and thus transform as such, 


i.e., if u = u(x) with inverse x = x(u), then 


W 


§' (u~b) 


Interaction Rate Densities 


Binary interaction description: 


s(x(u)-a] [J ® | 


Consider a beam of monovelocity particles perpendicularly incident on a 


"thin" 


slab of stationary targets. Let 


particles beam intensity (in em? 
target density (in en 4) 


target slab thickness (in cm) 


sec 


~L 


) 


Expect that the i-type particle-target interaction rate is proportional to 
NXI 
and define the “microscopic cross section" for an i-type binary interaction 


by 
Interaction rgte of 
type i per cm of thet = oN XI 
target slab 


Let Fy ~ i-type interaction rate density (in en sec), For the slab experi~ 


ment, 


aut ea: 
FE, =¥ (ON XT) No, I 


Cross section and flux: 
For the slab experiment, the "macroscopic cross section" is defined by 
ary = No, (in em”), Note that for a monovelocity particle beam, I = Nv. For 
the case of monoenergetic particles, define the particle flux by $= Nv (in em? 
1 


sec) and the interaction rate density expression is generalized to 


F,= &@ em ?sec™+ 


for a monoenergetic (but, not necessarily monodirectional) particle distribution 


in a field of stationary targets. 


Cross section energy-~dependence: 

Fy = No Nw for the slab experiment. Consider varying v and also varying 
N~ 2 such that Nv is v-independent. WN is certainly v-independent. Usually 

find that Fy = F,(v) implies o, = 9,(v), i.e., the "interaction area" varies 


with the particle approach speed. 


Example - neutron-nucleus resonance interaction 


oy, (y) 


The case of moving targets: 
Gonsider a particle with velocity v which interacts with a target with 


velocity V. 


Precollision conditions. 


In many cases, find that oF depends only on the "relative approach speed" 


eee 
9, = 6; (u- ¥) 


Consider the i-type interaction rate density of a particle distribution n(v) 


with a target distribution N(v). Define £,@W) by 


£,~@™) @yvaty = i-type interaction rate density of particles in dy 


at v with targets in ay at V. 


£,@W) avaey 


9, (lv - vidly - vinq@wd*w ary 
In these terms 
= 3.43 
F, * ff £,(v,W)d-vdv 
F, = Sf 0, yew!) |y-vin@ @arvaty 
Average macroscopic parameters: 
Average cross section and integral flux - 


It is often useful eefek to define the integral particle flux by 


d= unwary = Nv 


and, to then express the i-type interaction rate density by 


The averaging process for the cross section has thereby been designated as 


Sf oftv-w)|y - vinq@awa?vaFy 


J vay) dy 


Examples - 
1. Case of stationary targets, i.e., N(v) = N&(V) where N is the total 


target density (in em”), 


— Nf o,(w)vn@w)a?y 


f mw) ay 


vo, (v) 
ae aaa oe 


Vv 


Case of stationary targets and monoenergetic particles, i.e., 
n(v) = a, (8) S(v-v,) where n, (2 is the particle direction-of-—travel 
distribution. 


ae = No, (wv) 


Case of stationary targets and a classical equilibrium distribution 


of particles, i.e., 


3/ 


nw) = NGB)*/? exp{-mv?/2kr] 


< 
i) 
i 
~ 
~ 


= Cc a 
Z, = N= No 
i Vv 1 
Case of a + cross section and a classical equilibrium particle 


distribution. 


Z = NS. 
1 


Uncollided Particle Flux and Mean-Free-Path 

The discussions employed to define the microscopic and macroscopic cross 
sections G, and z,) can also be used as the basis of a first particle trans- 
port discussion. Specifically, consider the determination of the uncollided 
particle flux in a medium that is "illuminated" by a monovelocity beam. For 
simplicity, consider a beam perpendicularly incident on the plane surface of 


a homogeneous "half-space." As illustrated, the incident beam 


AN 


Ele) 


aon acer ener mn Oh I's 


intensity is I, the medium surface is at x = 0, and the uncollided component 
of the particle flux is denoted ® (x). The thin slab of medium between x and 
xtdx may be considered to be equivalent to the thin slab used in defining Oy. 
Confining attention to only the uncollided particles, the beam incident on 
this thin slab (at x) is o, (x), and, if there are first collisions (inter- 
actions) in the thin slab, then the change in the uncollided flux, d@ (x), 


across dx is given by 
dd, Gx) la ©, Getdx) - & (x) 
and is a negative quantity. The first-collision rate density at x is given by 


FG) = 2, (x) 


where the total cross section, Des is used because all types of first- 
interactions are counted. The first-collision rate per wnit area of the 
thin slab at x is F, (x) dx. First-collisions is the only mechanism for 


decreasing the uncollided flux, i.e., Fy (x)dx has the same magnitude as 


d® (x). Whence, 


do, (x) = F, (x) dx =- zr. © (x)dx 
d® (x) 
a a. o (x) and % (0) =I 


The solution of the uncollided flux relation is 


~2 x 
® (x) =TLe 


Mean~free-path: 

From the uncollided flux result a clear physical meaning of the 
total cross section can be deduced. The form of the result is immediately 
comparable to the radioactive nuclide decay relation, N(t) = Noexp[-atl]. Re- 
call the definition of the decay constant, A. Specifically, dt is the proba- 
bility that a radionucleus decay in the next time increment At (if At << 1). 
The statistical meaning of a. is, in analogy, given by: Bax is the probability 
that a particle experience an interaction (any type) in the next Ax of particle 
travel (if ze dx << 1). Moreover, the equivalent "as good as new postulate" 
is: For a particle traveling in a medium, the probability that it interact 
with the targets of the medium in the next Ax of travel is independent of the 
distance traveled before that interaction. Moreover, F, (x) = C3) only 
describes the statistical average of the first-interaction rate density at x. 
The statistical description is given by the binomial distribution b&(Z Ax). 
Note also that the description in location of first-collision, x, can easily 
be changed to a description in terms of the time (the particle spends in the 


medium), t, using the relation x = vt or Ax = v At. 
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Another route to deduce the physical meaning of a (which agrees with the 
above discussion) is to determine the average (or mean) distance (x) of par- 
ticle travel in the medium to the first-collision. Using the derived o,), 


the rate of first-collisions in dx at x (per unit area of this slab) is given 


by 


WE x 
By (x) dx = 2,9, (x) dx = pes e dx 


and these first-collisions are clearly generated by particles which have 
traveled a distance x in the medium before the first-collision. The mean 
value of the distance traveled in the medium to the first-collision is thus 


-=,x 
x u1 e ¢ 


-r,x 

t 
[ae e dx 
Oo 


The denominator yields the total rate of first-interactions in a unit area 
of the entire half-space and must be equal to I, the incident beam intensity, 
to have a steady state (t~-independent) problem. The numerator is equal to 


1/t,. Whence, 


and the immediate interpretation: ze is the inverse of the mean~free-path 
(of travel, x). The previously noted statistical meaning (viz., B Ax is 
the probability of interaction for a particle travel distance Ax) is certainly 
consistent with this mean-free-path interpretation: 1/2 Ax is the number of 
Ax distance increments required to "on the average" have one collision, which 
implies Ax/Z Ax, or 1/e,, particle distance of travel to a collision. 

It should be emphasized that in these discussions only the uncollided 
component of the particle flux has been considered. With reference to the 


illustration, 


eX ASS 


x 


if particles are scattered (or, reemitted in any type of interaction), then 
the particle flux at x, $(x), has contributions from those particles which 
have had one (2, (x)) or more (9,4), ®,(%) 5-4) collisions. 
The more complex situation of q = Go) is easily approached with a 

change of variable defined by 

x 

Q(x) = { E,Gx! dx! 
a 


and termed the “optical variable," or "optical distance." So-defined, the 
origin of 2 occurs at x = a. Often, the choice a = 0 is made and then 
&2= 0 and x = 0 coincide. Since the mean-free-path at location x’ is given 
by V/t., the integrand 2 Cx! )dx! represents the fraction of a mean-free- 
path across the slab of thickness dx'. Whence, %(x) represents the number 
of mean-free-paths between the plane at x and the plane at a (i.e., in the 
perpendicular, or x-direction). 

The differential equation 


do, (x) 


dn =- E.G) o, (x) and (0) =I 


is transformed to 


do (2) 


dn eatin $, (2) and (0) =T1I 


if a= 0. Note that 2, = Ox ()). The solution takes the form 


o() =1 Pag 


and, then transforming back to x yields 


9,(x) = 6 (Ux) = 1 eR) 


x 
= = ' ' 
, (x) I exp{ { t,x Jdx'"J 
° 
As an example of immediate application of this initial transport prob- 


lem solution, consider the uncollided particle flux which would emanate from 
a slab of thickness X which is uniformly and perpendicularly illuminated by 


a beam of intensity I. 


The slab is neither thin nor uniform, i.e., 
es = 0(1), or greater 


B= 2 (yy,z) 


The emerging, uncollided particle flux is given by 
o (Gy,2) = 1 eo Ky52) 


i.e., 


X 
@ (yz) = I exp[- | Ee Gey, 2)ax] 
fe} 


If the uncollided flux is the only particle flux component of interest, 
then the above result immediately yields the “shielding” properties of the 
X~thick slab. As an example of another aspect of the solution, consider the 
case where the slab composition is uniform and only the target density, N, 


varies, i.e., N = N(x,y,z) and 
E.Gsy.t) =o, N G,y,2) 


In this case, the result can be expressed as 


ete Dae ge le 
ae o, &, (Ky, 2) 


° 
and a measurement of o,@y.2) yields a determination of the x-direction, 
projected target densities in the X-thick slab. 

It is rare that only the uncollided flux is relevant. The multiply- 
scattered flux components either represent unwanted noise in the determina- 
tion of o> or they are of interest as being principally involved in the 
phenomena of interest. It is to the more complicated calculation of the 


entire flux, %, that we must eventually turn. 


ELEMENTARY TRANSPORT RELATIONS 


Particle Current 

Consider a differential plane area element d’s with unit normal vector @ 
indicating the spatial orientation. The elemental area is located at r where 
the particle density is n(r,v,t). In a small time interval At, the number of 


particles with velocities in d?yv which cross a’s is given by 


n(xr,v,t) ay Atysé as 


where veé = cos(v,6) with (v,@) the angle between velocity v and unit normal 
é. Note that a sign 


= dv 


convention has been thereby adopted relative to the arbitrarily chosen posi- 


tive direction of 6, i.e., 
v8 26 implies + flow of particles. 


The flow rate of particles in ay across d*s per unit area is given by 
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3 


lim ae n(r,v,t) dy At we a’s =n(r,v,t) dv wee 


Ato At ds 


With respect to the direction @, the positive (directed) component of 


the particle current, J t), is defined by 


and is equal to the rate of particle flow per unit area across a surface 
perpendicular to @ at location x at time t; counting only particles with a 
positive component of v along the direction é. 

Similarly, with respect to the direction @, the negative (directed) 


component of the particle current, J_(z,t), is defined by 


J(z;t)= = | veé n(r,v, t) dev 
ve8<0 


and is equal to the rate of particle flow per unit area across a surface per- 
pendicular to @ at r and t; counting only particles with a negative component 
of v along the direction é. 

Clearly, the net rate of particle flow per unit area across a surface 


penpendicular to @ at r and t is given by 


J,(t) - J_(@,t) = [ue n(z,x,t) dv 
Defining the particle vector current by 


Jt) = {z n(r,v,t) dv 


the net rate of particle flow per unit area across a surface perpendicular 
to an arbitrary direction @ at location rand time t is @-JS(xr,t). 

Case of an isotropic distribution: 

"Isotropic distribution” means that n(r,v,t) = n(x,v,t) when expressed 


in velocity Cartesian 3~space. If Cartesian coordinates are 


2~3 


used for the v-integration and v, is chosen in the é-direction, then 


stro too too 
J, (z,¢) = | | ( v n(zv,t) dv_.dv dv, 


a co 6} 
The Cartesian n(v) here actually depends on |v only. Thus, it is more 
convenient to employ the illustrated spherical polar coordinates, with @ 


the 


polar axis, for the integration, i.e., 


a/2 27 2 
J, (x,t) = ell | v cos@n(r,v,t)v” singdvdedg 
0 0 0 
which reduces to I, t) = = N(x, t) v(x,t). 


The Point Kernel 
Consider a point (at zr)» monoenergetic (speed Vy)s steady source of 
particles emitting one particle per unit time isotropically into an infinite, 


homogeneous medium. Let 


2-4 


n (2) = the steady uncollided par~ 


ottertan Sf ae. TOS ticle distribution resulting from the 
4 
, 
a b Ba point source (i.e., the steady state 
ep 
©) point kernel), 


the uncollided particle density =f ny (y) d>y, and 


Nj 


I,@ = the uncollided particle current =/ va (ey) dy. 


Intuitively, 
1/v 
N (x) = exp[-2_R] 
2 4nR be 
J, = A exp[-Z.R] @, 
4mR 
and 
Lv, 
n (r,¥) = exp[-2_R] 6 (v-v_ é,) 
° Ar? t oR 


As an example of the use of these results consider a configuration volume 
element ar in which the rate of particle scattering collisions is Eo(x)d°r. 
Presuming that the scattering is isotropic, the volume element dr appears 
as an isotropic source of strength Eonar. 


“a 


TL ae 


we Na 
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Denoting dd, (2) as the uncollided particle current at R due to particles 


emanating from @r having had a scattering in ar, 


aJ,(R) = B,0(x)a°r 


i exp[-E,R] @ 
4nR2 t R 


This result is now used to develop a heuristic derivation of Fick's diffusion 


law. 


Fick's Law of Diffusion 


e 8, 


1 


dt, = 2 _O(xr) 
P - 4tr 


Let dd, represent the particle 
flow density rate thru as due to 
particles which had last collision 


in dr, te. 


In terms of the illustrated spherical 
polar configuration coordinates (i.e., 
with the origin at a’s and polar axis, 


z, opposite to 6). 


exp[-E,r] cos6r” sinddrdedd 
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extent, 


With the (unrealistic) presumption that the homogeneous medium is of infinite 
1/2 27s 


oO 
J, = { { { Te (x) exp{-Z,r] cosOsingdrdedp 
0 60 0 


For algebraic simplicity, suppose that $(r) = $(z) and employ the Taylor series 


expansion 


= Oo. 
a(x) = @(0) + Gar (rcos8) +... 
to obtain 


rs co /2 36 ; 
J Eno: { { {o(0) + Gas (reos8) +. . .] exp[-E,r] cos@sin@drdé 
0 0 
With the (sometimes unrealistic) presumption that over the r-range (0, a few 


1/Z,) 00/dz is essentially constant and defining UL = cos6, 


z 1 
d= = c [ [o(0) + &, ru] exp[-Z 1] udrdu 
0 0 
= 8 ph i 3¢ 
Pe tg OR es “Nae 


a more general (and accurate) result is 


=i ~1ip@s 
34@) Zo) 55 Ox r 


where, Xg is the distance measured along the + é~direction, D is the diffusion 


coefficient 


ler is the "transport cross section" (is equal to a. for the case of isotropic 
scattering and weak absorption, qs << rS)> and r is now the location of the 
surface thru which particles are passing with reference to an arbitrary origin, 


i.e., 


Note that 


@ + Sr) = I) - IL) 


a6 
=P te. 
ae 


ff] 


and that the gradient of 6(r), denoted by Vé(x), is defined such that 


4 3 
@ + Vo(r) = ge». 
es 


Whence, the relation 
J(x) = - D V(r) 


which is Fick's law of diffusion. 


The operator V: 


—_ aa 
‘ 


iy) é 


Consider the orthogonal set of coordinates (x) 5X5 »%Xq) at the point P 
with the respective unit vectors (,,8,,8,). Presume that the x, have dis- 


tance dimension (e.g., cm). Let F(r) be a scalar function and note 


2-8 
p) 
1 2 
r) 
+ G >» X3 + o(x*) 
3 


But, (8, 58,585) is an orthonormal set and F(P') can be reexpressed 


a oF 4 oF 4 aF 2, , 2 
F(P') = FCP) + Fe, &, + Bx, e, + 3x, e,) x + 0(x*) 
where x = x8, + 8, + 48° Whence, 
F(P') ~ F(P) OF a OF «a OF na. , 
ries ee Sey e + Bx, °2 + Fe 83) ée + O(x) 
3 
and 
lim F(R')- F(P) _ OF _4, 
x20 a Dug =é VF (x) +X 5%q) 
where the operator V is defined by 
ae 220 in 3 a 3 
Vo =@ 4 +6, 7-+6 
pf ox, 2 ox, 3 8K, 
Examples -- 
1. Cartesian (x,y,z) 
Aa 
z &z 
-4 a 2 a 2 
Ve ta +t <, oe + ee 
x 
2. Cylindrical (r,9,z) 
A 
z ey os 
ce 
aoe 9 13 3 
; a ea. 2S ey ao 
: i ron . “9 5 36 * a, az 
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3. Spherical polar (r,6,o) 


é 
' rt 
' 
i res 
| 
16 m™ A 
mr ! ey VaQ we +e,4 2 46,4 p) 
1 i * &, Or or 36 > rsind 3¢ 
-6—-——-l-- -- 
re 


Continuity Relations 


Monoenergetic particles: 
Consider the balance in population of monoenergetic particles in volume 


V enclosed by surface 8. 


N(x,t) = particle density, 


J(z,t) = particle current, 


s(r,t) = particle source rate density. 


a { N(r, t) ar = rate of change of the particle population in V, and has the 
following contributing terms: 
1. sf s(r,t) der = rate of change due to particle sources. 
v 
ig f u(x,t) o(z, t) ar = rate of change due to particle absorption by 
Vv 
targets. 
3. —- f J(r,t) ° é d’s = rate of change due to net outflow of particles 
gr S cs gutriow 


through S. 


where a is the outward~directed unit normal vector to S. 


Note that by Gauss’ theorem, 


Thus, if the particles are assumed to be stable, 


= f NG, t) dr = - fV-+i3@,t) dr 
Vv Vv 


- f{ 3,(z,¢) 0(r,t) art f s(x,t) dr 
Vv V 


and, since V is entirely arbitrary 


NCL) 2 + FGyt) - ByCet)O(et) + s(z,t) 
which is the particle "continuity relation." For this assumed case of mono~ 


energetic particles 


ee ON(r,t) _1 30(r,t) 
; ot v ot 


which yields an expression for particle continuity in terms of the 2 dependent 
variables ® and J. 

Polyenergetic particles: 

In the more realistic case of particles with a range of kinetic energies, 
it is of interest to develop the particle population balance in E-space as 
well as in r-space. Here, considering only the E-integrated population balance 


in r-space, 


oN(r, t) 


sie me Ve SG,t) - 2, (x te(et) + 9(z,t) 


which is identical, in form, with the monoenergetic particle relation. 
Here, 

N(z,t) = f n(x,v,t)d-v 

(x,t) = fyn(z,y,t) dy 

It) = fya(z,y,t)d4y 


s(x,t) = f s(x,v,t)d>v 


and Da is the averaged particle absorption cross section as discussed in 


previous lectures. Moreover, in many cases, 


ON(r,t) __9 [O(r,t)}] _ 
att 7 


where Vv = f vn(x,y, t) av. 


tt _ 
n(x, t) 


In these lectures, then, the general form of the particle population continuity 


relation is expressed as 


ye Ot) = - Vs Jt) - 2 (x,t) o, t) + s(z,t). 


The diffusion approximation: 


Fick's law of diffusion has been derived in the form 
J(x,t) = - D(x, t)VO(x, t) 
This is an approximate relation between J and $¢ which is valid if: 
1. Distance to boundaries (or, to inhomogeneties) is much larger than 1/t.. 
2. There is negligible variation in 30/9 x, over distances of Lt, i.e., 
slowly varying $(r), or equivalently EA<<t,. 
3. There is negligible variation in $(t) over periods which are small 


relative to particle transit times. 


The "diffusion approximation" to the particle population continuity relation 


is thus 


4h 


= O(x,t) = V9 + D(z,t)vo(z,t) 
- £,(z,t)(z,t) + s(z,t) 
Note that, for the case of D(r,t) = D{t) the operator v7, defined by 


a 


is relevant, i.e., the diffusion equation assumes the form 


< |e 


= (x,t) = D(t) Wace, t) 


7 oo, (x,t) (x,t) + s(r,t) 


Examples -- 


1. 


Cartesian coordinates (x,y,z) 


9” 9? 
Pay a ee: 
ox oy oz 


Cylindrical coordinates (r,6,z) 


2 2 
2202 8. yh FL 3 
ie x or (r ax) * +2 38 + Je? 


Spherical polar coordinates (r,8,9) 


2 


se 8 oe 
Y= 2 ox (x 


3 
3 tz 
r 


Viscosity and Thermal Conductivity 


Viscosity: 
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Consider a fluid characterized by n(r) = n = constant, and v(x) = v(z) =v (ze 3 


From the Fick's law derivation 


= 


dd, = ee exp[-2.r] cosOsindrd6do 


Let P,.8, denote the net flow of particle momentum transfer through as per unit 


area. Then 


00 TT 
~ _Ss as _ J 
Pay a { { exp[ ir] v,(zcos8) cosOsinOdrdé 
00 
Note that it is assumed that © is approximately constant even though v, is vary- 


ing and that the $-integration averages ve at each value of @. 


25 om OO PTT es dv, 1 
Puz 2 { { exp[-Z.r]}v(0) - (reose) (Go), jcosSsinGdrdé 
0 0 
i.e., 
Vx 
Paz 7 neo 
1 z 
where L = ome md = Dm@é is the fluid "coefficient of viscosity." In three- 
t 


dimensions, the result is 


Thermal Conductivity: 


Consider a fluid characterized by n(r) = n = constant and v7 (r) = v"(2). 
Let Q denote the net rate of particle kinetic energy transfer through d’s per 


unit area, i.e., 


tm TT =z 
Q. = f | exp[-Z,r] vy’ (reos6)cos@sinbdrdé 


Zz 4 
0.0 
Moreover, if the particle distribution is locally in classical equilibrium, 


then 
PP 
3 mV (z) = 3 kT (z) 
and 
3 ie 
22 = i 
Q, Z KES [ | exp[ tr] T(rcos$)cos@sinédrdé 
= ~ St 
Q= -KD) 
L 2, 3 
where K = 7 Osa = > kw is the fluid "coefficient of thermal conductivity." 


t 
Qa, is identified as the "heat flow density." In three-dimensions 


Q= - KV 


MONOENERGETIC PARTICLE TRANSPORT RELATIONS 


Point Source Kernels 

To motivate a careful understanding and analysis of the particle trans~ 
port process, consider the details of a particular problem--Determine the 
particle flux, O(r), resulting from a unit strength, monoenergetic, steady 
point source of particles located at r' and emitting isotropically into an 
infinite, homogeneous medium. Two solution approaches have already been 
discussed in these notes. It is the relation between the two solution re- 


sults that is addressed here. 


‘“ flux cbservarrorns 
TF Pore 7 
ae f 


a“ 


0 woe / oo 


pe SOKrCE ae 
\ 


v\\ 


The uncollided particle flux, ®,@), has been previously discussed and 


is given by 


For reasons which will become apparent, denote this particle flux result 
by G(xr'), dee., 
G (xz, 2") = the uncollided particle flux at r due to a unit strength, 


2 monoenergetic, steady, isotropic point source at r' in an 
infinite, homogeneous medium. 


As a first detailed problem employing the derived particle diffusion 
equation, consider the unit strength, point source problem. Specifically-~ 
Find the particle flux, $(r), due to a source at r' using the steady state 
diffusion relation. In order to gain obvious symmetry, move the origin 
of spherical polar coordinates to the point source. Then @(r) > (r) and 


the diffusion equation can be stated 


id d 


SFE (Der) - o(r) = 0 


except at r = 0 where there is a suitable source condition. The "diffusion 


length," L, has been defined by the relation 


M 
ole 


The general solution is 
o(ey wah orb, od tlh 
r xr 
and the boundary (and source) conditions are 


(I) lim $(r) = 0 
mo 


(IZ) lim 4x? (-D =) #1 


ro 
Whence, the solution 
-r/L 
e@) -—— : 
4nd? 


Spacial moments: 


As expected, [2pm 4nx°dr = 1 which leads to the general spatial (r) 
moment relation 


<r> = \ o EO(r) 4nr¢dr = (nt1)11" 
0 
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@.g., <tr> = 2L and <r*> = 6L7, Note that <r7> # <r>?. 


Green's function: 
As for the uncollided flux problem, there are useful reasons for denot- 


ing the particle flux result by Gr"), i.e., 


Gy @ 5") = the particle flux at r due to a unit strength, monoenergetic, 
steady, isotropic point source at r' in an infinite, homo- 
geneous medium, where the diffusion approximation (equation) 
has been employed to determine the result. 

From the &(r) result obtained placing r' at the origin of r~coordinates, it 
readily follows that 
1 


Ear Teepe exp(-|r - r'l/L] 
4nZu* |x - "| 


Gy(z.zr') = 

which should be compared with the uncollided particle flux Green's function 
ty = 1 await 
G,(z.z') = —~——~> _ exp (-lx - x [/A 

an|z - x" 
where A denotes the mean-free~path, i.e., A= it... 
Use of Green's functions: 
The two Green's functions, G(r,r'), just described are certainly different 


in their functional form, i.e., 


~ 1 .-R/L 
G Re 
and 
~ 1 -R/A 
Sm RFE 


where R is the distance between source location, r', and flux observation 
location, r. However, both results are trying to describe the particle flux 
at r due to a unit strength, monoenergetic, steady, isotropic point source at 
x’ in an infinite, homogeneous medium. The relation between these two re- 


sults is considered later. Here, consider the general use of Green's 
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function results to determine particle flux due to complex source distribu- 
tions. 

Green's functions are generated by solution of linear analytical problems 
where the inhomogeneity (or, source) is singular - e.g., in the present prob- 
lem, point sources are considered, i.e, s(r) ~ S(r - r'). Since the problems 
are linear, the solution for a more complex source is given by the "sum" of 
the relevant Green's functions weighted with the distributed source density, 
i.e., the particle flux, %(r), due to a distributed source density, s(r'), is 


given by 
O(r) = [ ewe sears’ 


if G(r,r') is the Green's function which describes the physical process cor- 


rectly. If diffusion is the accurate transport description, then 
O(r) = [ sepx eG ar" 

(Or, if free-flight transport is the process, then 
O(x) = [ sez" yet! a%r! 


Relation between free-flight and diffusion: 

Using a method based on calculation of successive “orders-of-scattering" 
a relation between free-flight transport and the diffusion approximation will 
be developed. The purpose of this development at this point in these lectures 
is two-fold: 1. To qualitatively indicate the general idea of a transport 
calculation. 2. To show the inadequacy of knowing only the free~flight and 
diffusion description in studying general transport problems (and thus moti- 


vate a more careful consideration of the transport of particles in matter). 
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Calculation of the uncollided particle flux - 


If s(x') represents a particle source distribution, then 
= ' ' 35.1 
© { G (Ez) s(x')d’r 


is the uncollided particle flux generated by the source distribution. More- 


over, 


F(x) = 2,0) 


is the first-collision density rate. 
Calculation of the once-collided particle flux - 
If scattering is isotropic, then DoF, G@)/2, represents an isotropic 
source distribution of once-collided particles. Whence, 
z 
(= { G(x x!) z Fy (x')d3xrt = { G(x x") 26, (z')dex' 
is the once-collided particle flux. Moreover, 


Fy(z) = £,0, (x) 


is the second-collision density rate. 
Calculation of higher-order~collided particle flux - 
In general, 


oi = | G (erro (eae! 


is the m-collided particle flux (i.e., if scattering is isotropic and G, is 


particle-energy~independent). In such cases then, 


o@ = a { { ine [ eGeztdG, cet .2")+ 26,2 2D) (4D) 


a'r' d2x",. 3p) 
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Calculation of particle flux - 
If $(z) represents the particle flux (i.e., including all orders-of- 
collisions) due to the source distribution s(x"), then 
© 
a= J e&@ 
m=0 


There usually is some integer, M, such that 
M 
(1) * J ot) 
m=0 


and it is reasonable to presume that M = O(E,/Z,)- Two limiting examples 
have already been considered, i.e., 


1. Case of high absorption, i.e., 


Za >> zy implies M» 0 and 


M 
ney *m© * Mo * { Go (eoz")s(z')der! 


2. Case of weak absorption, i.e, 


a << a implies M>> 1 and 


M 
ne o> eo) * | Gy G@z')s(x')d’r! 


i.e., the diffusion approximation. 

If M is between these two limiting cases (or, if scattering is aniso- 
tropic, or, transport is particle-energy-dependent) then the orders-of- 
scattering technique is either tedious or impossible to apply and a direct 
consideration of the complicated particle transport process is required. 
The remainder of this section of these notes is devoted to development and 
solution of a specific form of the Boltzmann transport equation with appli- 


cation to relevant physical and engineering problems. 


Monoenergetic Particle Transport Equation 


Presuming that the particles of interest are essentially moving 
points (i.e., with no relevant external or internal structure), the distri- 
bution independent coordinates are location, r, velocity, v, and time, t. 
Alternatively, and more useful for the present discussion, the independent 
coordinates can be chosen as r, energy, E, direction-of-travel, Oy and t. 
In order to somewhat simplify the development, understanding, and applica- 
tion of transport relations, an often-used approach is to separately con- 
sider the independent variables (or groups of variables). In this section 
of these notes, the variable dependencies (x,8, t) are considered and thus 
monoenergetic particle transport is addressed. The next section treats 
the E-dependence and is devoted to particle E-spectrum determination as well 
as energy deposition (in matter) considerations. 

In the most complicated monoenergetic particle transport problem, 
n(r,%, t) is the relevant particle density distribution, i.e., 

n(x,0, t) = expected particle population per unit volume (e.g., per 

cubic centimeter) at xr, per unit solid angle (e.g., per 
steradian) at 2, at time t. 
The number of particles in volume V with direction of travel in the solid 
angle W at time t is given by 
| | n(r,Q,t)d?2 dy 


Vw 


where all the particles have the same energy (i.e., truly monoenergetic) 
or, for some reason particle energy is an irrelevant variable and thus par- 
ticles are counted independent of their energy. 

Consider the particle population change in an arbitrary volume V sur- 
rounded by a closed surface $ with local, outward—directed, unit normal e. 
Consider only those particles with velocity directions in the incremental 


solid angle d?2 about &. The rate of 


change of the expected population of relevant particles (i.e., those in V and 
479) is 
9 A 
d7Q Be | n(z,2,t)d?r 
Vv 


and has the contributing terms: 


1. dQ @,+viin(x,0, t)d?s = net rate of loss of particles from V with 
velocity directions in d° at time t due 
to flow through S. 


de (x, t)vn(r, a »t)d°r = rate of loss of particles in V with pre~ 
collision velocity direction in d°Q at time t 
due to binary collisions. 


v 
3. aa | d’r | e(x, te (x, tsQrmnye, (x, t)va(r, 9", t)d22" = rate of gai 
v 


4m of particles in V with velocity directions 
in 47Q at time t due to the production of 
secondary particles caused by collisions of 
particles in V with all possible precollision 
velocity directions. 
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directions in d?Q at time t due to independent 


4, d?#Q | s(r,0t)d3x = rate of gain of particles in V with velocity 
Vv 


particle emission sources. 


In term (3) two new functions describing the particle interaction process are 


introduced. They are: 


e(x, t) = 
f(x, 38190) = 


It should be noted 


the expected number of secondary particles released per 
collision of the particles of interest (e.g., c = m/e, if 
scattering is the only collision type from which particles 
emanate) 


the normalized secondary particle direction of emission (Q) 
distribution due to a precollision particle with velocity 
direction 2' (i.e., fd @ is the fraction of particles 
emitted into d?Q about 9) 


that term (3) includes the contribution of particle scatter- 


ing from all direction increments into d’Q about {% and all these scatterings 


are described by binary (point) collisions. With intent, the continuous devi- 


ation of a particle direction-of-travel which would result from long-range 


force fields has not been included in this initial, simplified discussion. It 


should also be realized that there are numerous situations where long-range 


forces are of negligible importance and the relations being developed have 


meaningful application to such cases. 


Gauss's vector field theorem applied to term (1) yields 


{ € -win(r, fi, t)a?s = | Veviin(r,G, t)d?x 


s 


Vv 


Moreover, including the fact that v and Q are xr-independent, 


| @-viin(z, 0, t)a?s = { vi-Va(r,8,t)d3r 


s 


v 


Combining the contributing terms (with proper sign according to gain or 


loss rate) and realizing that both d?2 and V are arbitrary gives the Boltzmann 


transport equation, 
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“sf a(x, Gt) = - vieVn(r,f,t) ~ ve, (x, t)n(z,f,t) 


+ c(z,t)v2, (x,t) | £(x, 032) n(r, 8c) d2a" + s(z,%,t) 
4n 


In terms of the particle flux, 6(r,2,t) = vn(r,2,t), the transport equation 


is 


dik 


gt 0(z,8,t) = - A-vo(e,A,t) - 2, (zt) oz, 8, t) 
+ e(x,t)Z, (x,t) | £(x,t38'>9)6(z,8',t)d2Q' + s(x,6,t) 


In this monoenergetic particle formulation of the relation describing the 
transport process, it makes little computational or conceptual difference in 
using n or % as the dependent variable. Both formulations will be employed - 
the choice depending on the nature of the problem. 

Conservation relations: 

Consider the operation { a2 applied to the particle transport equation. 


Term-by~term the results are 
[ @awa.e aa = { n(z,@,e@a = <2 wex,t) 
| vn (x0, t)a22 = ye [ incesl,eyata = VeS(x,t) 
[vt cen(e,d,eya'a = v2. (x, t)N(x,t) 
| cvE, | £(Q'90) nf" a2n"a2a 
ws [ evza@ [ 2@ fy atoarer = cvE,N(z,t) 
{ s(z,9,t)d?2 = s(r,t) 
Whence, 


se N(z,t) = -V-I(E,t) + [e(e,e) ~ LIVE (e,tNGt) + 2(e,0) 
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an equation which describes conservation of particle Xr-space population density. 
Note that multiplication of this conservation relation by the particle mass, m, 
yields an equation expressing conservation of mass (i.e., the simplest of the 
equations of fluid dynamics). 

More relations of fluid dynamics can be generated from performing various 
8-space moments of the transport equation. Such relations are far more mean- 
ingful when applied to the energy-dependent transport relation. However, the 
general idea can be briefly illustrated here. For example, consider the 
operation | a?Q ny, applied to the particle transport equation. The left~hand 


side of the relation is 


BG Ay Gdn Naas De 1 Be (a vvey ae 
{ W 5e n(xr,Q, t)d2Q mV 5 } Qa(x,Q, t)d?Q 


= av se 2 tN(E,t)] = go (a (cyt )N(E,E)] 


and represents the rate of change of the net particle momentum per unit volume 
of r-space. The right-hand side terms include discussions of complexity un- 
warranted at this point in these lectures. The resulting relation expresses 
conservation of particle momentum and is thus an “equation of motion" for the 
"fluid" of particles (i.e., a form of the Navier-Stokes equation of fluid dy- 
namics). 
Analytical Resolution of the Boltzmann Equation 

To illustrate techniques for attempting analytical resolution of the 
Boltzmann transport relation just derived, three approaches are outlined here. 
For algebraic "simplicity" attention is restricted to a special class of prob- 
lems. The Boltzmann equation we are considering has already included the 
simplifications of: 


1. No long-range forces and short-range forces described by binary 
interaction cross sections. 


2. Monoenergetic particles. 
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3. Cross sections independent of direction of particle travel @). 


4. Stable particles (i.e., no particle population change due to 
disintegration). 


To these add the following restrictions: 


Ds Homogeneous, steady state target medium (i.e., c and ze are constants 
and £Q! 38) is r and t~independent). 


6. Isgtropic emission of particles following a particle collision (i.e., 
£E(QQ'2Q) = 1/47). 


7. Particles in steady state (i.e., n = n(x,2), or ® = $(r,9)). 

8. Plane symmetry. 
In detail, the restriction problems with plane symmetry means that (x,2) > (x,8), 
where the position coordinate, x, and the direction-of—travel coordinate, 9, 


are defined by the illustration 


A 
2 


é 


Clearly, in order to actually achieve plane symmetry all medium boundaries or 
discontinuities would need to be in the form of infinite planes perpendicular 
to the x-axis, and particle sources would be required which are in the form of 
x-planes with emission-direction symmetrically distributed about the x-axis. 

The restriction of plane symmetry as well as the other seven simplifica- 
tions listed above are never completely descriptive of a real problem. How~ 
ever, techniques for analytical resolution of more general transport problems 
are usually based on the ideas outlined here for this greatly simplified 


case, 
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For this case, the dependent variable appears always as vwn(x,@) and 
therefore the flux (9) equation requires less symbols to express. In addi- 
tion 

§-V0(x,8) > cos6 52 4x8) 
T 
| £(G'98) (x, 0")d?a" + $ [ #¢x,0")sin8"a0" 
4 ie) 


is more naturally expressed in terms of the direction-cosine, Ll = cos@, i.e., 


@ + o(%,y) and 
B16 +n <2 oxy) 
ox ou 


1 1 
| £od7Q" + = [Seu ay" 
4v -1 
In these terms, the relevant Boltzmann equation is 
5 chy +1 
Hg OGw) + E0Cow =e f eeu" aut + sGaw) 
-1 
It is to this relation that the three solution approaches discussed in these 


notes are applied. 


Approach 1. Eigenvalues and functions of the Boltzmann operator: 
Consider the homogeneous form of the Boltzmann equation (s = 0) and in- 


clude effects of s by boundary conditions. Thus, 


3 ch 1 
Gi ay t 2.) OCxu) = 7 O(x,n')du' 
1 


Translational invariance of this equation suggests solutions of the form 
O(xu) = ¥(2,u)e*/* 


where ¥(£,u.) must then satisfy 


er +1 
(en ues EL) ¥@.u) = > | ¥C2su" du! 


od. 
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It is found that acceptable values of 2 fall into a discrete and a con- 


tinuous spectrum. The continuum in & is for real & in the range 


= vt, SRE W/E, 


The discrete values of %, are 2 = + L with c < 1 implying that L is real and 
c > 1 implying that L is imaginary. 
The eigenfunction set {¥(2,u)} is orthogonal in the sense 
+L 
| p¥(2,w¥C2', udu = 0 for 2 # 2! 


<1 
which can be expressed as 


' 
(CxSner> 2 and 2' discrete 


1 
fora ¥(2' ,y)dy = 
se |e,ect-2), & and %' continuous. 


The general solution of the homogeneous equation is in the form 


L/E 
t 
2(xu) = AGL)YGL, we + act ¥CLwe™™ + ih A(R)Y (Rue *! Mag 
2 Ee 
where the set of expansion coefficients {A(2)} is evaluated by boundary condi~ 
tions with the aid of the orthogonality relations. 
Approach 2. Expansion in direction-dependent functions: 


Assume that the set of functions {2, Gu} is complete for -l <y<tl 


and use the expansion 
O(x,y) = } Fy (x), (u) 
in the homogeneous form of the transport equation, i.e., 
d ck 7 1 
enti? oy 2 ' 1 
Tuy) Ge FLO + EM GWF, OO] = > LF g(x) My Cu" )du 
-1 


A particular example of a useful function set is the Legendre polynomials 
{Pp(u)} which, in addition to being complete, are orthogonal on -l Su <+1 


in the sense 
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+1 2 
| PoGud Py 1 Gudu = Teri See 


and obey the recurrence relation 
(24L)Po 41 (H) ~ (241 )uPy (u) + Py, (u) = 0. 


For simplicity of physical interpretation, redefine the expansion coefficients, 


Fy (x), such that 


Hu) =] 272 6, (x)Pp (u) 


and note that, with this redefinition, 


+L 
99 (x) = 2m { Gu) P,Gdu = fO(x,u)P,(u)a?e 
-1 


For example, 


®, (x) fox, u)d?Q = (x), the particle flux 
® (x) = fo(x,wud?2 = T(x), the particle net current 


The transport equation becomes 


yi 2241 


d ck “1. 
tan (UP, Gu) ax Oo (x) + EP (H) 8) Gx) = —z7 Oy (x) | Po (u' du") = 0 


#1 


which, using orthogonality and the recurrence relation, takes the form 


& 
YC Py Hq FAP, GD gate + 20 + 1)P,QDE,,60 


ee cL, (x) 690] = 0 
or equivalently, 


d d 
, CL Fe Opin) + GH) Fe Op, Ge) + (QUILL) (x) - cE, 0, (x)5)0] Py) = 0 
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Recognizing that the set {P, (u)} is linearly independent, 


a (x) + (241) gE 0), G0) + (2RFL)E,D) (H) ~ cB, 8 (by, = 0 


= 9 
dx “2-1 R41 
i.e., an infinite set of coupled differential equations in terms of the ex- 
pansion coefficient functions, %, @)- 

The Pyvapproximation: 

A method for truncating the above set of differential equations is to 


presume that (x) = 0 for 2 > N and only retain the first N+ 1 equations. 


As an example, consider the P,-approximation equations 


do, (x) 
ae + (-e) 2,6 (x) = 0 
dd, (x) 
a + 32.) (x) = 0 
which imply 
2 
Ae Bee ite Gai tS ae 
3, ax* i 


Defining the diffusion coefficient D = 1/3E, 5 and identifying (x) and ®, (x) 


as particle flux and current respectively, the above relations are 


2 
D aye - (1c) E(x) = 0, the diffusion equation 


I(x) = - p 42) 


is ae 23 Fick's law of diffusion 


Note that the solution of the homogeneous diffusion equation is of the form 


(x) = B ent ty yp ota 
where the diffusion length, Lh is defined by 


41/2 
D 
a feos] 
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Approach 3. Moments method: 


Including an extraneous source, s(x,y), the transport equation is 


(ia “= + 2) OG, = ar om ie $(x,u" du? + s(x, y) 
1 


Using the Legnedre polynomial expansions 


ow = J 2 a cor. 
g 

s(x) = J) a 89 (x)Po(u) 
g 


reduces the transport description to the set of differential equations 


L 4 4 R41 _d 
Weel Se Pa) FTO Tae ee) + Gedy) 2,90 Cx) = sy (x) 


Define the particle flux distribution moments and particle source moments by 


i dux™ Pe (u)9(x,u)a72 


-o 00 


amd n, 
Pon = { x O) (x)dx 


400 foo 
Bane te | x89 (x) dx { dxx” JP, Qu)s(x,n)d?2 


For example, 


boo too 
mn 
2=0: oe | dxx” f6(x,y)d?2 = | x(x) dx 
00 on 
Pom m, : ‘ ¥ A 
Therefore . " <x >, i.e., the configuration moments of the particle dis~ 
fers) 


tribution. 


m= 0: $= {- dx f Pp Quo(x,u)d2Q = ("aw (x) 0 (x) dx 


co) oO oO 

Therefore e = <P Cu) > » i.e., the direction-of-travel moments of the particle 
00 

distribution. 


In order to obtain an equation relating the {o)} with the set {s,,} per- 


form the operation 


00 
{ x ™ (equation above in o 


oo 


Pea? Spar Spex 
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and note that 


‘e do al 
& mm = m, to _ { m1 = 
im “i * dx = (x 0) eo nN mx , (x) dx By, ml 


Whence, the relation 


[26 i+ (241) 6 


A-e5y ES pq * Spm * Feet (MQ gl, m1) 

In many problems, the source moment set {so} is given, or can be obtained 
from a given source distribution s(x,y), and the particle flux solution is 
sought. The results, here, indicate a systematic, algebraic technique of ob- 
taining as many members of the particle flux moment set {o,} as is practical, 
or desired. From the particle flux moments either certain physical descrip- 
tions are immediately evident (such as problem scales), or an attempt can be 
made to artificially construct the particle distribution, 0(x,u). It should 
be carefully noted that although there is an infinite sequence of algebraic 
moment relations, they are coupled such that % | depends only on 4, m1 and 
Pou, mL’ Thus, a systematic progression of obtaining moments given {3 is 
generated. The set {4} depends only on {so} and is thus presumed given. 


The following matrix indicates the systematic progression in the calculation 


of moments. 


ke e) ZL 2 eave &-1 L SL 
mR 
a ? ® oo 1, P60 Poo direction 
moments (presumed 
: %, given) 
2 9a, 
eae J ae Peel, med 
m o, —w 
= ® 
mt &n 
. configuration 


moments 
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Application and Validity of Analytical Methods 

In terms of the analytical approaches just described, attention can 
now be focused on application to physical problems. Such considerations 
help to clarify the details of the methods as well as illustrate their limi- 
tations. 

Use and validity of Py approximations will be illustrated by applica~ 
tion to the problem of determining the reflection properties of a homogeneous 
half-space. Use of the moments method will be illustrated by appiatation to 
the problem of determining the effects of source emission distributions. 

Reflection properties of a half-space: 

Consider the problem of a radiation flux illuminating the plane surface 


of a homogeneous half-space. With reference to the illustration, a 


% 
A 7 
Vaccum a Hlomogene ous wrediunn 
pon be od 
oO 


homogeneous medium in the positive x half-space is illuminated by a particle 


flux with a prescribed angular description, i.e., 


6(0,u) =¥@) forOg us +1 


The emerging radiation flux, $(0,u) for -l £0, completely describes the 
reflection properties of the half-space as dictated by the parameters of the 
medium and the incident illumination. Note that we have already discussed a 
more restricted descriptive quantity, R, the reflection coefficient, defined 


by 
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Certainly Rcan be evaluated via radiation flux illumination and emerging 


descriptions, viz. 
=-1 0 
J,(0) = 2 | u¥(u)du and J.(0)+-25 LOCO, py) du 
0 -1 
Thus, the problem is fully directed toward determining the emerging radia- 
tion flux, 6(0,y) for -l<u<0. 


The Pyvapproximation is summarized by the truncation 


N 
zs 2ntl 
O(x,u) = > an 8) PL) 
n=0 
and the equation set (for isotropic emission) 
n4 @ (x) + (ott) & 0 (x + (2nH)E (ed) (x) = 0 
dx n-l dx n+l t no’ n 


with n = 0, 1, 2,..., Nand % (x) = 0. 


N+1 


Pj -approximation~— 


The particle flux is given by 
=—t 3 
GW = 7 Oo) +7 Gey 


and the relevant equations are 


$= 4) + T.U-e) Gx) = 0 
4 g(x) +35. O (x) = 0 
dx %0 aes Baca 


We have previously discussed the manipulation of these relations to generate 
Fick's law and the diffusion equation. Let us now approach the solution 
of the problem: in a manner analogous to finding eigenvalues and eigenfunctions 


of the Boltzmann equation (Approach 1, page 3-13). Specifically: Translational 


invariance suggests the form 


eG) = ¥ exlh 


Whence, 
RE (i-c) ¥)(%) - ¥#(® =0 
-¥ 7 
Go) + 3k ¥ 0 


i.e., two linear, algebraic, homogeneous equations in the two unknown ¥9(2) 


and ¥, (2). The condition for non-trivial solution (where the trivial solution 


is ¥y (2) = ¥, = 0) is 
ad (1-c) -1 


-1 322 


which reduces to 2 = Ly where 


= [3(1-c)]74/2 a 


The eigenvalues *L are the P,~approximation to the discrete pair of eigen- 
values tL discussed on page 3-14. For a c<l problem, the eigenvalue spectrum 
for the exact solutions (Approach 1) can be compared to the eigenvalue spec- 


trum of the P,~approximation (Approach 2) via the illustration 


Alppro ach / Approac 42 
Lack Sold Z - Approx a atior 


The P| -approximation eigenfunctions are easily extracted from the above 


relations. Specifically, 


¥ = 1 


% (2) 322, 


yields ¥, GL) /% GL £1/3L,2,- Whence, acceptable solutions of the 


4) = 


exponential form are 


0 1 
O(x,u) = {1 +75 Ele 
an Liz, 
and 
¥y(-L,) +x/L 
O(a) =S ef-s eule + 
lt 


i.e., a general solution of the form 


x/L tx/L, 


ir 


1 S 1 
$(x,n) = A[L +7 yle + c{l - syle 
are Lee 


In the half-space reflection problem, the boundary conditions are 


(1) lim (x,y) = 0 


X Fo 


(II) A condition descriptive of the entrance 
particle flux ¥(u). 


Boundary condition (I) implies C = 0 and thus 


1 -x/L 


L 
®(x,u) = A[l +> ule 
or ha 


In particular 


2(0,n) = AIL + ul 
ive 
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Clearly, the application of boundary condition (II) cannot accurately describe 
an arbitrary illumination ¥(y), O< u<+1. To satisfy the P| ~approximation, 


the linear y-form is required and only the very restricted case of 


¥u) = all ++ uw] forogusetl 
tA, 


would give direction-of-travel pointwise agreement. It is therefore usually 
required to develop less rigorous boundary conditions which are, at least, 
intuitively satisfactory. An example is to establish entrance particle current 


agreement, i.e., apply boundary condition (II) as 


+L a8 
20 { u6(0,u)du = 27 a uUY(W)du = 
0 
which is merely a statement that J,,(0) =I. Using this arbitrarily chosen 


condition yields 

ree: Rae 
wll + 2/3L, 1 
and the problem of characterizing the emerging radiation is solved, i.e., 


- I 
a(0,y) = WI + 275, 37 {1 + W/L, 2] 


for ~1 <u <¢+1 which includes the emerging direction contribution. In pass-~ 


ing, note that 


1/2 


= [3G-c)] 


and thus the expected result that half-space reflection properties are only 
a function of ¢ and, for example, are independent of Zee Rewriting the emis-— 


sion distribution as 


00,u) = A[1 + BG-c)u] for -1 <¢uso 
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illustrates another Limitation of the P, approximation. If 1/3(i-c) < 1, i.e., 
if 0< é@< 2/3, then $(0,y) < 0 in the range of emerging directions which is 
a physical impossibility. The case of ¢ = 0.5 is illustrated to indicate the 


clear limitations of the P,-approximation applied to this problem. 


ue 


PO) 


ee a 


“a 


a 


2 


ae Maumtrnatron 


| 

| 

\ 

. 
Emerging se 
we | 


of 


Finally, with perseptive disbelief for cases where ¢ < 2/3, the calcula~ 


tion of the reflection coefficient proceeds via 
0 
| I 
w{l + 2/3L,2,] 
-1 
to the familiar "diffusion approximation" result 


J_(O) = -2n HL + u/L,2, du 


- zr - 
J_(0) L- 2/3L,4, 1 2D/L, 


Ro Same eee eS Oe 
I 1+ 2/3L,2, 1+ 2D/L, 


P,-approximation—- 


The particle flux is given by 


ps 
w= O@ +e {out 2 om BY-H 


and the relevant equations are 


GE a) + 2,G-c) a(x) = 0 


ad 
dx dx 


d = 
2 = ® (x) + Sz. (x) = 0 


Same methods as used in the P| ~approximation give now 


2E,(1-c) ¥g(2) - ¥,@) +0 = 0 


- ¥y (2) + 322, 4 (2) - 2 ¥o (2) 


0 - 24%) + 52E .¥, (2) 
where 
a(x) = ¥ (aye /® for n= 0, 1, 2 


For non-trivial solution, 


£2, ce) -1 (0) 
-1 SRE. 2 =0 
0 -2 SkE 


t.e., 
~poes 
15(1-é)2 Da - [5 + 4(1-c) ]22 = 0 
which reduces to 2 = 0 and +L, where 


2 


4 1/2 


hn mi 
WL" Gaatpl %, 


Oye) + gad (x) + 32, (x) =0 
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The eigenvalues +L, are the P,~approximat ion to the discrete pair of 
eigenvalues +L. Comparing the discrete eigenvalues of the various approaches 


for ac <1 problem yields the illustration 


O B-approx. Dor 
Ba Ze Bpperow. 
O freacxH i 
| 
—Lz the 
Bet —__#- Gre 
-£ ~L, ae <= tL, EZ 


The P,~approximation eigenfunctions are extracted from the above re- 


lations and yield the two ratios 


= AZ (1-e) 
3 dod l 
- 2 P22 (1-0) -$ 


‘ie ie! = 2 2 ~ . 
Whence #y GL.) /¥5 (425) = 4150, c) and ¥, G13) /¥5 (45) = 3L5 2,1 @)/2-1/2; 
acceptable solutions for the particle flux take the form 


: _ 3,2 , 342 1, be 
Ox) = aha f1 + 3L, 2, (i-e) + 5345 2 (1-2) ~ rE) GS ue ple 


In summary, the P,-approximation solution to this problem yields a better 
asymptotic eigenvalue than does the P,~approximation (i.e., Ly is closer 


to L than is L Moreover, a more complete angular distribution is 


a" 


generated. However, because the additional eigenvalue, 2 = 0, is of no 


use in satisfying boundary values, there is no more flexibility found for 
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more completely specifying the illuminating particle flux ¥(u). For this 
reason, the next approximation used in problems such as the one considered 


is the P,-approximation. 


3 


P,~approximation~- 


The particle flux is given by 
3 2 1 
O(x,u) = ad O00) + oF a ou + z ®, GEG UW = = 


at Pare 
+ Ba@e wv - dy 


and the relevant equations are 


sf &@ +4-¢) &@w =0 


dx 

<2 ox) 42-4 8 (x) +35. = 0 
dx 0 dx *2 tts 

2 fo(x) +358 2 6,(x) + 52, $5(x) = 0 


d 
3 a 4 Gx) +7 ze , (x) = 0 


x/ 2 


Using & (x) = Ye forn = 0, 1, 2, 3 yields the set of relations 


&E. = c) ¥, (2) = YO + 0 + 0 = 0 


~Yy (2) + SREY C2) - 24, (9 + 0 = 0 
0 = 24 (2) + SREY, CX) - 3¥(9 = 
0 + 0 - 3%, (9 + 7 RE, ¥3 0%) = 0 


and the resulting condition for non-trivial solution, 


RE (1 - ¢) -1 0 0 
al 32 -2 0 
t 
= @ 
0 -2 5 oe ~3 
0 (0) -3 7 2%, 
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d.e., 
105(1-c)2"n0 é [55(1-c)+35]2722 +9=0 
Whence, 
2.2 _ 354¢55(1-c) 3780(i-c) 4/2 
whe = “Bioti-<) ie) ive 2 
[35+55(1-c)] 
Note that 


[35455(1-c)]7 = 1225 + 3850(1~c) + 3025(1-c)2 


and defining e(c) by 


3780(1-c) 


e(c) = 7] 
[35+55(1-c) ] 


it is easily concluded that, for 0<¢ <1, 0<e(c) <1. Moreover, asc +l, 
é(c) + 0. To indicate further relevant numerical details, consider a case 


with c sufficiently close to unity such that the approximation 


fr - e(ey 2 2 a - hee) 
can be employed. Thus, 
age = SAG) [ys(d-e(c)/2)] 
Whence, the solutions 
L = tL, where L3 = [ say + + ie ye 


353455 (1-c) t 


- 1/2 scl 
3 3 


2 = +L3 where LZ = [ 


The eigenvalues +L, and +L7 are the P,~approximation to the discrete pair 


3 3 
of eigenvalues +L and the continuum of eigenvalues between ~1/k, and +1/ Bes 


respectively. Comparing the eigenvalues of the various approaches for a 
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ce <1 problem now yields the illustration 


O F- approx. H A- plane 
A 4 approx 
QO Fy7 appre. 
o Exact 
z 
Z , n 2 
ie ee Ae | +L3 te +L, +h; 


The P,-approximation eigenfunctions are extracted from the above relations 


and yield the three ratios 


¥, (2) 
CS il eae 

¥, (2) 

a 2 ea 2 

Way 7 PE Be) - 5 

¥,(2) 

Ec ar Oe ORE 
Ty. 2 ge Be oe ae 


Acceptable solutions for the particle flux take the form 


Y, (+L.) Y, (4L.) 
= ENe3 2535 Bed: 
O(x,y) =A +3 ¥oCt5) wuts ¥Cz,) Gu 3 


se ¥,G2,) é 2 3 i. as 
¥oGL,) z 2 
W, (+12) ¥, (+13) 
Ps i233 2 53" 32, 
+ A°/1 +3 “+5 pq Ge -pP - 
b GL) ¥yCh3) 2 2 


¥, (412) ~x/L5, 
B35. 3: 8 fag) 
ae ¥9 C23) Ge -3 5) . 
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Exclusion of the otherwise acceptable solution forms which vary as exp(tx/L,) 
and exp (+x/L3) has already invoked the boundary condition that $(x,y) > 0 as 

x +o, The two coefficients, A and A”, are determined solely by conditions 
descriptive of the entrance particle flux ¥(u). Thus, like the P,-approxima- 
tion, the present solution yields a progressively better asymptotic eigen- 
value (i.e., L, is closer to L than is either Ly or Lo)> and a more complete 
angular distribution is generated. In contrast to the P,-approximation, the 
additional eigenvalue pair 2 = +L and their associated eigenfunctions provide 
additional flexibility in more completely characterizing the illuminating 
flux. 


In the P,~approximation $(0,u) has the form 


3 


o 2 3 
6(0,u) =A (Cy + Cu +C, uo + Cc, uw) 


sas pps & ey S53) 
+A (cg + cru + C5 ywo+ C3 uo) 


where the c. and ome are known functions of the half-space medium parameters. 
The cubic p-form is required and only the very restricted case of 
YQ) = A, + A + A? + AGW for 0 <p <tl 
o TAH + AS 3H gus 


and the fortuity to find that there exist an A and A” which satisfy the 


simultaneous over-specified set of relations 


- 


AC, + A‘C 


i> 
i] 


0 0 19) 
AY = AC, +& cy 
Ay = Ac, +4 cy 


Ag = Ac, +X C3 
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yields direction-of-travel pointwise agreement between the solution $(0,u) and 
the illumination boundary condition y(y). It is, therefore, again required 

to develop less rigorous, but intuitively satisfactory, boundary conditions. 
Using the same approach as was applied to this problem under the P,-approximation, 


an example is to require 


+1 +1 
20 { wo(O,uW)du = 27 | uY(WdH = I 
0 0 
and 
+1 +1 
21 | y°9(0, udu = 27 { wy du = # 
0 


which is a statement that not only J, (0) = I, but that the average value of 
u-cubed is also in agreement with the illumination conditions (i.e., a one- 
step better description of the actual illumination). 


Using these two arbitrarily chosen boundary conditions yields 


L B°L-o°H 
T 


a a8” = a8 
aoa db aH - @1 
T of -~a°B 
where 
a = ¢+4c¢ 446,44 
0 1°22 "5 %3 
os c+ ocr¢icre2cz 
Gre ee ka a 
1 2 2 
8 Fea» Sia i er ae 
- h Laer iy oe ee RR a 
8 Tote t+FCe+ Hc 


Thus, the problem of characterizing the emerging radiation is solved. In 


addition, the calculation of the reflection coefficient proceeds via 
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0 
J_(0) = - 27 | H®(0,u) du 
“1 
which yields 
B°L-O°H | — OH - BI 
J_(0) = % oe ap + OB as 
where 
~ 2 1 2 
Cos tz, - FC, 
a, Po eee ane: ae ey ae 
Se RG. a ee eg g/g Se 


The reflection coefficient is now easily expressed in terms of half-space 
material parameters (the a,a°,B,B, @ and a’) and illumination description 
(the I and H) as 


J_(0) 


2 ~% BI- aH | O° on - Br 
id I ag* ~ a8 I ap ~ a7B 


An alternative algebraic reduction-- 
Beginning with the Py approximation set of couple differential equations 


for the functions { (x), n = 0,1,2,...,N}, i.e., 
n4-6 (m) + (nt) © 0 (x) + (Qn + DE (1-c8_)9 (x) = 0 
dx “n-1 dx ntl t no’ n 


with n = 0,1,2,...,N and Poa = 0, another, more direct, algebraic reduction 
is possible for problems such as the one just detailed. The approach follows 
the general path of that covered in Approach 1 (page 3-13). Specifically, 
translational invariance of the equation set suggests solutions of the form 


—x/2 
o,(%) = ¥ (Le 
In fact, using this form of solution actually generates a recurrence relation 


for the "eigenfunction" set {¥ (2), viz. 


n ¥-1 ~ (2n +1) (1 - 6) REY C2) + (at1)¥ = 0 
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As written, the relation among the YO is reminiscent of the equivalent 
recurrence relation for the Legendre polynomials. Moreover, since these 
solutions will eventually be linearly combined (summed) with coefficients 
chosen to match arbitrary boundary conditions, an arbitrary normalization, 
say ¥ QQ) = 1, may be employed. The desired function set is then generated 


by the recurrence relation, viz. 


¥ (2) 


W 
te 


¥,Q) = (i-e) 22, 


y 


3 252 
¥, (2) gle) a - 


Nie 


_5 3,3 _ 54+4(1-c) 
¥,(2) = FG-c)e"ne - AEN a2 


35-4_n\p4p4 _ 35455(1-c) 252 | 3 
¥,2) = Fc) REP oe MEE +s 


etc. 


and, as expected, bears a resemblence to the Legendre polynomial set. The 
Py approximation is accomplished by using only n = 0,1,2...,N and setting Vey (2) 
= 0. 

It is no surprise that the set {¥ ()} so~generated is identical to the 
functions found in the earlier discussion. What, perhaps, is a surptise is 
that the previously determined P,-approximation secular equations (condition 


N 


for non-trivial solution) is precisely the same as the relation Ved 


In fact, comparing the secular conditions and the present requirement yields: 


(2) = 0. 


Pj~approx. secular cond. = [2 Yo (2) = 0] 
Po~approx. secular cond. = [6 ¥3(%) = 0] 


P- approx. secular cond. = [24 ¥, 0 = 0] 
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It can be shown, with some algebraic effort, that these results generalize to 


Py approx. secular cond. = [(N+1)! Fen (O = 0] 


Since it is relatively easy to generate the set {¥(D}, this is a useful 
relation. The existence of the relation lends a sense of consistency to this 
method of truncating the expansion of 9(x,1). 

Source emission distribution effects: 

Consider the problem of a single, plane radiation source in an infinite, 
homogeneous medium. We shall employ the moments method (Approach 3) to esti~ 
mate, ina relatively simple manner, the influence of the source angular 
emission distribution on the resulting radiation flux. 

The flux and source moments have been previously defined as 


o = | xP (U) Ox, u)dxd72 


nm 
-o 4 


San ft x (u) (x, W)dxd?2 


“oo AT 
If the plane source is positioned at the origin of x, then s(x,u) = s(u)6(x) 


and the set of source moments is generated by 


1 
Sim * 27845 : P,(@sQdu = sd 
“1 


For the purpose of detailing a specific problem, consider the case where s(u) 


is a symmetric, quadratic function of the direction variable H, d.e., 


sty) =At cy? for -l Susg+i1 


In terms of the Legendre polynomial expansion of the source emission distribu- 


tion, 


1 5 3.42 L 
sw = Gr So tae 82 Gu - DD 
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Clearly, 
at 5 
Be ee Rg. 5 Bo) 
aids, 
C= or 82 


With the (angular integrated) strength of the source, Soe fixed, the limits 
on sy which generate a positive-definite (s(W) 20 for -l <u <+1) source 


distribution are found by requiring A > 0 and s(+1) > 0, viz. 
A20 + s, ¢ 2s 


s(#1) > 


Vv 
° 
+ 
a 


Whence, “Sp /5 £8 s 2s,/5. The limiting cases, depicting the range of 
symmetric, quadratic source emission distributions of physical relevance are 


illustrated. 
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With 89 and 85 (within the required limits) chosen, the source moment matrix 


is established, i.e., 


ae 
0 L 2 3 
nm 
, 0 2 0 5 0 
l 0 0 0 0 
2 0 0 0 0 
3 0 0 0 0 


The radiation flux moment set ia Pe is generated by the previously derived 
> 


set of relations 


m 
~ QntL ee Ree - (tL) Oa eal * Ee -c6) 8, - Sam 


The resulting low index corner of the flux moment matrix is 


n> 
0 1 2 hes 
; 0 “0 0 22 
2, (i-e) z. 
1 50 
1 0 5 Gag +2 ay) 0 
3 
= 
2 89 
2 Gs 2 85) 0 
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One technique of using the determined flux moments for practical calcula- 
tion purposes is to formulate a guess of the functional form of the radiation 
flux and then apply the exact moments in the "optimization" of that guess. 

For example, in the present problem spatial symmetries imply that (xy) = 
$(-x,-1). Moreover, previous lectures suggest an exponential form for the 


x-variable and a polynomial form for the p-variable. Whence, the guess 


O(x,y) = A[l + BP, (yu) + cp, (ule! # 


where the upper and lower signs are applied to x > 0 and x < 0, respectively. 
The choice of using four descriptive parameters (A,B,C,2%) is based on the 
intention of employing the four flux moments found previously (i.e., Poo? 


9 0 and wD for parameter optimization. In detail, 


+00 +1 
®. = 27 i { O(x,U)dxdu = 8TAd 
aL: 


7 ~, 


teo pt] 2 3 
® | = 27 { { x O(x,U)dxdu = 167AL 


02 
—_O =] 
Ho ehh 8 
®0 = 20 | | P, (1) (x,y) dxdu Be ALC 
=o 71 
to tL 8 2 
a, { { 2P (u) S(x,y) dxdy = 57 AL"B 
to Ay 
which yield the relations 
3 1/2 
wh ®o2 1/2 . L 00 
2809 An 299 


Using the flux moments derived previously for this problem in terms of the 


source emission parameters S85 and 85> 


ee 2 
oo Ga 
3a l-c Sp 


s 2s, )-1/2 
ca Xe of | 
l=-c 


R= 


81 I-c 


a 
it 
uw 
on 
pay 
{ 
ie} 
VY 


These results used in the guessed form for the radiation flux, 0(x,u), give 


a solution which is optimized in the sense that the moments ®oq° %yos %o and 


41 are exactly correct. 


A comment about diffusion theory-- 


If, in the above discussion of application of the moments method, S_ = 


0, 


then the results are descriptive of the radiation flux generated by a single, 


isotropic plane source located at x = 0. For this case 


ees Pree Ot 
‘00 uc) 02 3231-0)” 
from which 
® 
ees - 54 - 2+ =212 
00 32) Gc) 


where Ly represents (as previously) the asymptotic (actually, only) eigenvalue 


obtained in the P| -approximation. 


It is interesting to note that if a P| -approximation to the present 


problem is used, the result 
O(x) = 5 (x) =A exp[-|x|/L] 


is obtained. Whence, if <x>, represents the second x-moment under P,-theory, 


‘ / 
-x/L. 
2 [ xe Lay 2 
<x“%>_ = — = 2L 
1 ~a/Ly 1: 
e dx 
0 
Thus, <x2> = <x?> + That is, the rather surprising result that however well, 


1 
or badly, the P,-approximation yields valid representations of flux solutions, 


or parameters based on flux solutions, it always gives the exact second 
x-moment. There are certain physically interesting problems which can be 
demonstrated to depend primarily on <x?> (e.g., the probability that radiation 
will emerge from a finite slab having been generated in the interior) and, 
these problems are well-described by diffusion theory even when material par- 
ameters dictate that the P)-approximation is invalid. 
Numerical Resolution of the Boltzmann Equation 

The previous two sections of these lecture notes are devoted to the 
description and use of techniques which are primarily analytic. By "analytic" 
is meant that results may be formulated in general functional forms which 
indicate the influence of relevant parameters on the trend of solutions. 
Numerical evaluation of the analytic results is only required to visualize 
the details of these results. In this section, solution techniques are 
analytically generated which require extensive numerical computation to yield 
any meaningful information. These approaches are termed "numerical resolutions" 
and two of them are outlined here. Again, for algebraic and notational sim- 
plicity, attention is restricted to the special class of radiation transport 


problems which are described by the list of simplifications used in considering 
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analytical approaches (viz., items numbered 1 through 8 on pages 3-11 and 3-12 
of these lecture notes). The relevant form of the Boltzmann equation is (from 


page 3-13) 


os +L 
i ae OGyu) + BOG, = = O(x,u)du’ + s(x, W 
-1 


It is to this relation that the two numerical solution approaches discussed in 
these notes are applied. 

Approach 1. The discrete ordinates method: 

The primary objective of the discrete ordinates approach is to reduce the 
integral term of the Boltzmann equation to algebraic terms. The emission of 
particles in collision (e.g., scattering) is described only in the discrete 
set of directions-of-travel fu, n= 1,2,...,N}. With this procedure employed 
to reduce the integral term, it makes sense to use the same idea on the other 


terms in the equation. Thus, the transport equation takes the form 


N 
uD S 
Uy Ge On) + E,oGou,) = WO Gay -) + sQru) 
n=] 


for n = 1,2,3,...,N 


where the w, are "weight" functions suitably chosen with reference to the dis- 
crete directions HU, So as to optimize the summation approximation to the actual 
integral term. The set fu,» ¥sn=l,2,...,N} is termed the quadrature set and 
the choice thereof is important in developing a valid approximation to the 
solution of the transport equation. 

An example of the generation of weight functions is illustrated by the 
use of trapezoidal quadrature. In this case, it is presumed that $(x,u) varies 
linearly with j; between any two successive quadrature directions, say Hy 


1 
and Uy» i.e, 
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nol E 
= co) + ¥ 
som MMe Gl uyYa-1 a “a 2 


With this presumption, the contribution to the integral term, [oG,u*) au’, from the 


i a to is 
interval Ho-1 Hy 


ps S(x,n7)au? = + 
J * 2 
Ho-1 


Similarily, 


L 
pty P Cou) + F Gy yo Gu, 4) 
ntl 


ayace od 1 
SGuu du” =F Gy Heo) + F Gy) oGeu,) 
Mn 


and, these two u-interval contributions contain the only terms in which (x54) 


appear. Thus, the contribution to fOCx,u7)du7 from the quadrature direction My 
is 

F ggg Hy 7) 8G) 

2 S"ntL “n-1 ai 


and the weight function is identified as 


ws mu ) 
n 2°'ntl "n-1 


This result only applies to the case of trapezoidal quadrature and is based 


on the presumption that $(x,H) is as illustrated. Other assumptions 


B(x) 
© Fr u,) 


ae / LU Sa LS3 O 
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on the functional form, or moments, of ®(x,y) lead to different quadrature sets. 
It is instructive to examine solution of the homogeneous (s = 0) form of 

the discrete ordinates transport equation along the same lines as applied to 

the Py approximation. The equation is invariant with respect to translation 


of x which suggests the solution form 


oss ~x/2 
$(xu) = ¥@,u ye 
Whence, 
N 
CRE, 
ULE Yu) + 2, Wy -W(2su,-) = 0 


for n = 1,2,3,...,N 


which is a set of N linear algebraic, homogeneous equations in the N unknowns 


{¥@su.), n= V5.2, woe NI. 


For the purpose of clarifying use of this result, consider details of 


the problem characterized by N = 2, uy = Uy» and Wy) = Wy = 1. The discrete 


direction symmetry chosen is as illustrated 


where a = - 85, u, = cos8 = cos8,. The equation set now has only the 


a? Hg 
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two members 


CRE, 
[uy ~ 22, C-c/2) ]yCasuy) + = Ys) = 0 


u 
OQ 


cky 
[uy = 2 (1-c/2) 1¥ (2,5) a ¥(A,U) = 


Using H, = - Uy and requiring non-trivial solutions for FSB) and ¥(RUy) 


generates a secular equation with solutions 


RL=+b where L=+ 


The eigenvalues (+L) can be set at a desired value by suitable choice of the 
quadrature direction Uys For example, if Hy is chosen to be 1/73 (which 
corresponds to 85 = 54.7°), then L= L,> i.e., the P,~approximation result. 
In a later section of these notes the solution just obtained (i.e., for 
N= 2, Wy = Uys Wy = Wy = 1) will be applied to a specific physical problem. 
Here, it is of further general clarification of the discrete ordinates idea 
to consider the discrete u-dependence of the radiation flux for positive x, 


i.e., 


/L 


il 


O(x,u,) = A ¥CEL, ue 


/L 


OX; Ly) A W+L, we * 


A description of interest is the flux "angular spectrum" which can be expressed 


by the ratio (x, 1))/(x,1,), i.e., 


(x, H,) ¥rL, u) c/2 


8G) EGU) ~ Cee/2) = re 


which is surprisingly independent of the choice of equadrature (i.e., Vy) 


However, as expected, the illustrated graph of the angular spectrum is generally 
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greater than unity, approaching unity as c+ 1 and approaching unbounded 
value as c > 0 (which would yield no radiation in the H,-direction). 

To this point in this discussion, the discrete ordinates method has been 
presented as based on transformation of the direction-of~travel variable, u, 
to discrete form leaving the position variable, x, in continuous form. The 
motivation for such a presentation is primarily for comparison with previously 
discussed analytical resolution approaches. We now proceed to the equation 
forms which are usually numerically resolved with the aid of digital computers, 
i.e., a discrete treatment of both L and x. 

Discrete ordinates method with discrete position variable~- 

The radiation flux, evaluated at discrete values of both direction-of- 
travel and position, takes on the indexed form 2x su) which is reminiscent of 
the radiation flux moments em As is the motivation in the development of 
the moments method, the full discrétization of the flux specification is 
performed with the purpose of developing a strictly algebraic set of relations 
which can then be resolved by digital computation procedures. 

The operation applied to the x-variable which must be transformed to a 
discrete form appears in the "flow term," uy dé/dx. In discrete terms, a 


possible formulation is 
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49,1) Gry ot) - OGc,5 4) 
eT 


dx Xu 7 %n 


and the discrete ordinates equation can be expressed in the form 


O(x 1,54) - Ou) 
ml? “n mn = 
et a Ox HL) 


N 
ck. _ _ 
ae Fi > wi Ou) + sGr ul) 
n7sL 
for n = 1,2,3,...,N3 m= 1,2,3,...,M-1 


The notation Ox ou) indicates that it is appropriate to perform an averag- 


ing process over the x-interval xn to x in expressing such terms. For 


motL 
example, presuming that (x, 1) is approximately linear in x over the inter- 


val Xn me 3 would indicate 


< 
= aH 
nu) = PO ety) + OG HL) 
Hy 2 
It is unnecessary to make any assumptions regarding 8, that is the averaging 
process for radiation source terms since, s(x,U) is presumed to be a given 


function and the evaluation of s(x, su) should proceed via 


Xmt1 
= | s(x, uy) dx 
S(x,u,) = —2—-—8 
tl ~ *n 


The discretized transport equation has been reduced to the form of N(M-~1) 


linear, coupled, algebraic equations in the NM unknowns 
{oGq ou.) 5 n=1,2,...,N; m= 1,2,...,M} 


For each discrete direction-of-travel, Lae the radiation flux is illustrated 
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ANN 


e B46, fbn) 


As indicated, the positions x) and Xy are usually associated with actual 
physical boundaries (they are always associated with boundary conditions). 
The incoming radiation flux at the boundaries generate boundary conditions 


for the discrete flux of the form 
O(4) ou.) =¥,0q), n=No +1,...,N 
® (you) = ¥ C5)» n= 1,2,...,N7 


where the YG) and ¥ 052 are either given incident flux conditions or are 
derivable from fluxes other than those in the above boundary sets, and N° < N 
is determined by the largest value of n such that aH <0 (i.e., H, through 

through UH, are positive). With the N boundary 


u,,7 are negative and U 


N N°+1 N 
conditions specified there is a total of NM equations in the NM unknown 
(x Ju) and solution should be possible. 

A solution strategy for the discretized equation-- 

Since the flux induced emission (including scattering) and radiation 
sources are employed in the same manner in the iterative solution strategy 


outlined here, use the notation 
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which should not be confused with the same notation used earlier in these 
notes for the radiation flux moments (page 3-17). The discrétized transport 


equation, rewritten in these terms, is 


n,mtl nm ot - 
Ma ax, bs 2 CO, antl + ee Qaim 


for n =1,2,...,N; m= 1,2,...,M-1 


where ax = Xe 7 Xp and the boundary conditions become 
eT Yi Gy), n=N/2 +1,...,N 
aM = Y, Gy)» n= 1,2,...,N/2 


where for obvious symmetry advantages the quadrature set fu} is chosen such 
that N is an even number and Wy, = My» Uy = = Uyeperses By so =~ My /24i' Two 


forms of the discretized equation are useful, viz. 


1- 2 Ax /2u 


am 


® Se ee 
na, ortl 1+ 2 Ax /2u, om u/2 + u/Ax, 

‘ 7 L+ 2 Ax /2u, : . Qam 

am L- 2 Ax /2h n,mtL ri /2 - ui /Ax,, 


To describe the particular solution strategy presented here, it is useful 
to detail a problem which is, at least, specified by a given class of boundary 


conditions. To this end, consider the problem of a slab of uniform material 
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illuminated with a given radiation flux on the face x = 0 and occupying the 


region 0 $x § X. At the face x = X, a specific reentrant (or, non-reentrant) 


flux condition is imposed. The general problem and specification of position 


Homogeneous ned iu 


|: slab ee 
Eryven i Cven 
Jociden? teen Tran? 
Slew he axe Ax cond) P1244 

pee — afc ~/< <0 

Osasts | : 

1 1 | ERDESAL SOIT Reena ease _- non : ae 

Kp Es ig Re Xugep 
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mesh points (with uniform spacing Ax = X/(M-1)) is illustrated. With the 
previously mentioned condition of N even and a symmetric U-gradrature set, 


the boundary values at x = 0 can be expressed as 
= i n = N/2+1,...,N 


where (Ys n = N/2+1,...,N} describes the illumination radiation flux. The 
boundary condition at x = X is described when it is required in this formulation 
of a solution strategy. 

The starting point (or, zeroth-order interation) of the solution is to 
select (guess) a set {a3 n= 1,2,...,N; m= 1,2,...,M-1}. To signify the 
iteration index, denote this set {Q,,(0)}. Actually, it is presumed that the 
Som contribution to Oot is given and thus the selection is for the remaining 
radiation emission term (e.g., a possible guess is Q m6) = Samm + K). The 


iteration strategy is accomplished as follows: 


3-49 


1. Use the "forward progressing" form of the discretized transport 


equation, i.e., 


. . 1- 2 Ax/2u F n Q m6 
Tl, otL 1+ EAx/2u nm 12 + u/Ax 


for n = N/2 +1,...,N; m= 1,2,...,M-1 
Denote the results of this calculation om) signifying the results of the 
first iteration (note, only for the forward directions yo > 0). 


2. The last calculation of step 1 determine 


@ Ql), n= N/2t1,...,N 


which are the emerging (transmitted) radiation fluxes of the first iteration, 
forward calculation. The boundary conditions at x = X are applied to these 
flux values. For example, if an operator % describes the reflection proper~ 


ties of the space in the region x > X, then 
® Gh) = Olen DH, n = 1,2,...,N/2 


which provides an “effective illumination" on the slab face at x = X for 


continuation of the first iteration calculation. 


3. Use the "backward progressing" form of the discretized transport 


equation, i.e., 


ee 
P 1+ Boe /2U : r Qn (9) 
= T-.t&o2u E72 - ux 
nm 1 7 Ax/2y, nymtl z/2 H/Ax 


for n = 1,2,...,N/2; m= M-1, M-2,...,1 


Denote the results of this calculation ome which completes the first 


iteration flux results by supplying the components n = 1,2...,N/2 (i.e., 
HL < 0) at each position x The values of ® yA) for n = 1,2,...,N/2 are 


determined in step 2 and the Qa) have already been selected (guessed). 


4. The second iteration begins with determining an updated value for 


the set (Qa? i.e., with obvious notation 


Q dss + wo oy Wye (Oe ae D +O. DI] 


Steps 1 through 3 are then applied to generate the set {9 62) 3 WF Legescag Ny 
m= 1,2,...,M}, and, the iterative strategy is established. 

Illustrative examples of the "albedo" operator, a, are helpful in under-~ 
standing the role of boundary conditions at x = X in the iterative process. 
Three relevant cases are: 

A. A non-reentrant surface at x = X, i.e., the condition that eM =0 
for n = 1,2,...,N/2 independent of the iteration. For this case a is just the 
zero operator. 

B. "Specular" reflection at x = X is described by requiring that the 
emerging flux OM for (N/2 + 1) $n /N only acts a source for the symmetric 
guadrature direction Uyeion* Thus, independent of the iteration, the set of 


relations 


= for n = 1,2,...,N/2 
oM een *yet—n,M nee, , 


with each oO < 1, describes the reflection at the surface x = X. 
c. “Diffuse” reflection at x = X is described by requiring that the 
emerging flux oom for (N/2+1) ¢ 1 § N acts as a source of equal intensity 


for all the reentrant directions (i.e., Wa = 1,2,...,N/2). Thus, inde~ 


pendent of the iteration, the set of relations 
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N/2 
- £ = 1,2,...,N/2 
om 7 > ke PE eee / 
n-=1 


with each a, < 1, describes the reflection at the surface x = X. 
The solution strategy progression of calculation is illustrated for the 


three boundary condition cases at x = X for the particular choice N = 4 and 


M=5 ; 
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Approach 2. The Monte Carlo method: 


The method outlined in these.notes (pages 3-4 through 3-6) under the title 


"orders-of~scattering" suggests another numerical approach to solution of the 
transport equation. Recall that the successive scattering method is based on 
initially calculating the uncollided flux, o), by a convolution of the 

free-flight Green's function, Gy(E), with the source density, s(r*). This 


result provides a means of calculating the first-scattering density, i.e., 


L69)@), which, in turn, provides a source for the once~collided radiation flux, 


$, (x), which is calculated by a convolution of Gy (E) with Loo (x). The 
process is continued to determine the twice-collided flux, (x), the thrice- 


collided flux, $, (x), ete. The total radiation flux, $(r), is given by 
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where the value of M for valid truncation of the sum depends on the value of 
c. This method of determining $(r) is based on progressively following the 
expected trajectory histories of particles as they "migrate' through a medium 
from their point of "origin" at a source. The Monte Carlo numerical approach 
is based on a similar tracing of the path-histories of particles, but, in 
contrast to the orders-of-scattering method, the point-of-view is to simulate 
the statistical interaction sequence by a numerical analog experiment. 

Suppose that we have the means to generate random numbers, &, in the range 
0 <&€¢ 1. Another way to express this is that the probability of choosing & 


in the interval d& about € is independent of the value of &, i.e., if 


P(€)d—& = the probability of choosing the value 
of the number € in the range § to § + 
dg, 


then, P(E) = 1 for 0 <&<s 1. In order to use our random number generator for 
the purpose of simulating radiation travel and interaction, relations must be 
developed which couple the statistical aspects of the physical process to the 
probability distribution P(€). As an example, consider the free-flight travel 
of a particle in a given direction in a homogeneous medium from a point which 


we arbitrarily designate as the origin of distance coordinate r. If 


p(r)dr = the probability that a particle will suffer 
the first collision (after leaving r = 0) in 
the distance interval r to r + dr, 


then, it has been often used in these notes, that p(r)dr = ae exp(-2 r)dr. We 
would like to find the transformation r(E), i.e. & > r, such that P(E) gen- 
erates p(r). Thus, 


P(E) = p(r) =1 


dr 
ag 


for Og Es landr>0 
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Whence, 
ge] = ptr) = exp(-Z,r) 
E(@) = ~ exp(-E.r) 
Eic(é) = - &E 


Similar to the discussions on page 1-8 of these notes, a graph of rer(&) il- 
lustrates the idea of using a random number generator to simulate a physical 


process (in this case, free-flight transport). Note that random 


van 


—t 
Oo. ne oe 6 1.0 
numbers in the range 0 < § < 0.1 yield values of free-flight distance, r, in 
the range -£n 0.1 = 2.3 < her < ©; random numbers in 0.1 <& < 0.2 generate 
1.6 < yt < 2.3; random numbers in 0.2 < & < 0.3 generate 1.2 < ur < 1.6; 
etc. Since P(€&) is a constant, the first-collision density, which is pro- 
portional to the uncollided particle flux,is maximum at r = 0 and monotonically 


decreases with increasing r, viz. the decreasing exponential function exp(-Z ir). 


Expanding these ideas to describe particle emission direction probability 


densities, consider the case of isotropic emission, i.e., 
PCH, >) dudd = = dudd = Ea du = do 
. 4a 2 20 


where pp = cos@ is the polar scattering angle direction cosine (-l <u < +1), 
@ is the azimuthal scattering angle (0 < $¢ < 27), and identification of 1/2 
with du and 1/27 with d? is justified by the part of the normalization which 


comes from each variable, i.e., 


p(an = Fay 


it 


a(o)do aq ad 
P(uied) = plyq(o) 


We require the transformations y(&), i.e., E+pu, and o(n), i.e., n+ >, which 
yield the direction-of-emission distributions, p(y) and q($), from the random 


numbers € and yn. Note that 


PCE,n) = pCu, >) [sce | =1 


is simplified to 


P(&) = p(y) ie =p 
ain) = 4(¢) {g oa 


since ) and $ are independent descriptors of particle emission. Whence, 


Iss = pd -5 
ew) = F@4D 
ue) = 2E = 2 
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and, 


les = q(¢) -+¢ 


ne) = 5-6 


(n) 


it 
3 


a a 


In terms of the illustrated coordinates (with bases vectors ey, eos &,) which 


ra a 
relate pre-collision direction-of-travel, 2°, with emission direction, 2, 
Q. = sin®cosh 
Q, = sinOsingd 


2, = cosé 


and the direction cosines of a particle emitted from a collision are established. 
The values of @ = arccosy and 9 are generated from the respective random vari- 
ables € and n via the transformations u(E) and O(n). It should be noted that 
for the case of isotropic emission (and, only for that case) the coordinate 
bases set (8),8),83) may be chosen with a fixed relation to location coordinates, 


x. For example, using the illustrated conventional relation between 


A 
ez 
Q . = sin®cosd 
x 
Q. = sinSsind 
y 
2) = cos8 
z 
A+ 
Ke) 


Cartesian location coordinates (x,y,z) and spherical polar direction coordinates 
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(6,6), the emission direction, a, has the indicated components. 

Point source emission directional dependence, s(u,), also yields trans- 
formations between random variables and direction cosines. For sources, the 
coordinate bases set is always chosen with a fixed relation to location coor- 
dinates, r, e.g., 0858) as indicated in the illustration for isotropic 
emission. 

At each collision a decision must be made as to what particular inter- 
action type occurred. Suppose that there are I interaction types possible 
indexed by i = 1,2,..., I, and described by the cross sections zy. The 


probability that a collision (which has occurred) is of the i-type is 


Py = Z,/%, where te = > ort 


i=l 
The transformation of random variable, &, to interaction type, i, is relatively 
easily accomplished. With the definition of random variable limits 
i 
&(i) -> P 
j=l 
the random variable range 0< € < €(1) indicates that a 1-type interaction 


occurred; the range €(1) s € § &(2) indicates a 2-type interaction;...and the 


range €(I-1) ¢ € $ 1 indicates a I-type interaction. 


Application of Numerical Methods 


It is significantly more tedius to detail application of numerical resolu- 
tion techniques than is the case for strictly analytical approaches. Not 
only in the demonstration of validity, but also in just obtaining anything 
but trivial results there is usually a requirement of extensive calculation. 
This can be reasonably accomplished only on digital computers. 

Use of the discrete ordinates approach will be focused on the calculation 


of a homogeneous, half-space reflection coefficient and the procedures and 
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results may be compared with the Py approximation application to the same 
problem. The Monte Carlo approach will be discussed as applied to the deter- 
mination of the radiation flux generated by a point, isotropic source in an 
infinite, homogeneous medium. 

Reflection coefficient of a half-space: 

Consider, once more, the problem of a radiation flux illuminating the 


plane surface of a homogeneous half-space with a prescribed angular descrip- 


tion, 


$(O,u) = ¥Q@i) for O< p< +1 


The specific example of the discrete ordinates method discussed in these notes 
on pages 3~42 through 3-44 (i.e., N = 2, Hy =- Uy» Wy) = Wy = 1) can be 
applied to the determination of the reflection coefficient. Using discretiza~ 


tion of only the U~variable, the radiation flux for positive x is given by 


=x/L 
OOH.) = A ¥CLu,)e 
®(x,Uy) = A VrL ude ™!t 
where 
u 
L=+ none dea 
Zz vl-c 


and, the discrete angular-dependence ratio 9(x,Uy)/PCx,H,) is given by 


¥4L,uy) F 
e/2 
¥aGu) 
(1l-c/2) - Yl-c 


As is the case of the P, ~approximation approach to this problem (pages 
3-20 through 3-24) we usually face the situation of not being able to exactly 
meet the {illumination boundary condition with the required general solution 


form obtained for the approximation. For the Pj ~approximation it required 


the fortuity of a special linear U-dependence to exaxtly meet incident boundary 


values. Here, only the "monodirectional" illumination case of 
¥(H) = C S(u-u,) 


yields point-wise agreement between incident flux and required solution form. 
The alternative, practical approach of employing intuitive relations between 


the given ¥(u) and the required form for (0,4) can, again, be applied, e.g., 
+1 
J, (0) = 2m, 9(0, uy) = 20 { uY¥(u)du = 1 
0 


That is, the equivalent of the current condition suggested for the P,-approxi- 
mation. 
Using the arbitrary normalization ¥CtL Hy) = 1, the present general solu- 
tion form, and the incident current condition, 
2TH, Ac/2 


J,(0) oa came ernie Tae 
(1-c/2) - Vi-e 


and, the radiation flux for x > 0 is given by 


(l-c/2) - Y1-c i -x/L 
— e 


Tey, 


O(x, wy) a 
2 


J -x/L 


(x, Uy) = Tu, Le 


Moreover, the reflected radiation current at the medium surface, J (0), is 


then given by 


3_(0) = - 2m, 40,4) =Gas = Mire 


where the condition nn) has been applied. The reflection coefficient 


is thus 


70) Gea) - vice 
I c/2 


It should be noted that this result for R is independent of the particular 
incident flux boundary condition chosen to determine the radiation flux. 
This is analogous to the same situation found in the P| ~approximation (and, 


the P,-approximation) where only one useful eigenfunction allowed little 


2 
flexibility in solution form. It is, perhaps, surprising to find here that 
R is independent of the choice of Uy (or, equivalently, L). 

The illustrated curves of R as a function of c for the Pj ~approximation 
and the present discrete ordinates approach provide a comparison between 
the two methods for this type of problem. The complete breakdown of the 


P approximation for low vlaues of c is expected and is briefly predicted in 


the discussion on page 3-24. 


——— Pi/screr 
| ordinates, VaR 
aa a a Zo aperoxcina (oe | 
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Radiation flux from a point source: 

Demonstration of a complete Monte Carlo approach to solution of a given 
physical problem requires extensive numerical calculation. Such is considered 
outside the scope of these lecture notes. To indicate the steps required in 
an application, an outline of the sequential procedure is presented here. 
For the sake of notational simplicity, the simple problem of determining the 
radiation flux resulting from a point, isotropic source in a homogeneous, 
infinite medium characterized by the two interactions, capture (nonproductive 
absorption) and isotropic scattering,is presented. 

The sequence of decisions required to perform a Monte Carlo simulation 
of the point source problem is illustrated by the indication of a possible 
trajectory emanating from the source location which is chosen to be at the 
origin of r-coordinates, i.e., point 0. With reference to the illustration, 


the Monte Carlo approach steps follow: 


Two random numbers are generated, say Bo and Noe which establish 

the emission direction, Qos via the transformations derived on pages 
3-54 and 3-55. Specifically, Hp = z Eg7h and = 2m™M9- Note that 
the fineness of the mesh describing emission directions, Qos depends 
on the number of significant digits retained in the specification 

of the random numbers, e.g., random numbers generated as 0.00 to 
1.00 yields 101 possibilities for Uy and 101 possibilities for $5» 
whereas 0.000 to 1.000 yields 1,001 choices for Hy and $: 

One random number is generated, say Eo» which establishes the first 
free-flight distance-of~travel, Ry? via the transformation derived on 
page 3-53. Specifically, ERo =- & Eo Note that the location 

of the first collision, Eo» is easily determined by the relations 

Xp = Rysin&,cosd,, Yo = Rysin®)sing), and Zp = Rycos 8) (where oo = 
arccos Uy). 

One random number is generated, say Xg> which establishes the type 
of the first interaction. Using the notation discussed on page 3-56, 
let x(1) represent the ratio r/t, (rounded off to the number of 
significant digits in the random number generator). Then, following 
the ideas expressed on page 3-56, if % & X(1) the collision is a 
scattering and we proceed to step 4; if % > X(1) the collision is a 
capture and the simulated history of the first emitted particle is 
terminated (indicating a return to step 1 and emission of the next 
particle). In either case, the collision is recorded in a manner 
that both r 


—O 


lated particle trajectories have been completed, the sum of all the 


and ay can be recalled. Note that after all the simu-~ 


collisions which occurred in a volume d7x6 about Xo 


direction-of-travel in A7n about a, will establish the radiation 


with pre-collision 


flux at (£9 9M) via, Ee 9%p) is proportional to the collision 
sum. 

4. Two random numbers are generated, say Bl and Nye which establishes 
the first-scattering direction, ay» via the transformations uy = 
25 - 1 and o, = 2m, - 

5. One random number is generated, say ee which establishes the second 
free-flight distance-of-travel, Ry» via the transformation 2 Ry = 
- &nE, The location of the second collision, Ey: is determined by 
the relations x =X + R, sin8, cos, , ¥y = ot R, sin&, sind, , and 
2, = 2 + R,cos®, (where 8, = arccos Hl, ). Note that the simplicity 
of these expressions is a consequence of the presumption of isotropic 
scattering. 

6. One random number is generated, say Xe which establishes the type 
of the second interaction. Specifically, if Xy S$ X(1) the collision 
is a scattering and we proceed to step 7; if X 2 X(1) the collision 
is a capture and the simulated history of the particle is terminated 
(indicating a return to step 1 and emission of the next particle). 
Again, in either case, the collision is recorded such that ry and 
a can be recalled. 

7. Steps 4 through 6 are repeated requiring random numbers & and N 
to establish &, random number Ey to establish Ry; and random number 
Xz to establish the interaction type for the third collision. The 
location of the third collision, Zy> is determined by the relations 
X= Xy + R,sin®,cos>,, Yo * + Ry sin®,sind,, and Zo = 2y + 
R,cos@,. 

Thus, the Monte Carlo simulation proceeds following each emitted particle 


through its "life history" which terminates when the particle experiences a 


capture collision. As indicated previously, the radiation flux is determined 
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(after a large number of particle histories have been generated) by the rela- 
tion 
2,0(2, 4) A72A"9 is proportional to the sum of recorded 

collisions which occurred in rr about r 

with pre-collision direction-of-travel in 

area about 2, 
Generalization of the ideas to more realistic problems is straight-forward 
but usually algebraically complicated and tedious. Since many particle 
histories are required to generate useful information, the computational 
costs (use of a digital computer) is usually high. Several strategies are 


employed to emphasize the radiation of particular interest in a given problem. 


PARTICLE ENERGY SPECTRUM RELATIONS 


Classical Equilibrium Distributions 


Consideration of the E~dependence of particle distributions can proceed 
directly from the transport equation and such will be accomplished later in 
this section. Energy spectrum formulation is often more difficult to visual- 
ize than configuration spatial distributions-~thus, let us begin this dis- 
cussion with some of the more directly addressed problems. One of these is 
the derivation of thermodynamic equilibrium particle distributions. For the 
moment define "equilibrium for a particle distribution as: 

1) Steady state, i.e., t~independence 

2) Uniformity, i.e., r-independence 

3) Isotropic motion of particles, i.e., &-independence. 

Whence, with this definition, n(r,v,t) > n(v). 

Maxwell-Boltzmann distribution: 

The method employed here is based on a variational technique used in 
Maxwell's original derivation. Use Cartesian v~coordinates, i.e., v = 

= 92 gn on 


WysVy¥,) and let dn = av, Wx + ayy dv, * ov, dv 


Assuming that n(v) is separable in the form n(v) = n_(v_)n_(v_)n_(v_) yields 
xx’ yoy Z°z 


at Ee gh re gh Ci eh ge 
d 
n ny dv, ny qv, y a, dv, 


Consider a special change in v such that v is a constant, i.e., 


2 2 2 2 
vo tv +vo =v" = con 
‘- y Zz stant 


For such a change in v, 
1. dn = 0 since n = n(v) 


2. 2vdv_ + 2v_dv_ + 2v dv. = 2vdv = 0 
x x yoy ZZ 


Note that only two of the dv, are arbitrary (i.e., independent). With this 


change in v, 


L dn, 1 dn. 1 dn, 
aa dv, + aa 7a dv. + a. cage dv, = 0 and 
x x y y y Zz Zz 


2ov_dv_ + 2av. dv + 20y dv. = 0 for an arbitrary constant o. Whence, 
x x y sy ae 


1 da. Ll dn 1 dn, 
(i +" + 2av_)dv +(=— — + 2 Dav + (— = + 20v_)dv_ = 0 
n, dv. x Xk n_ dv. y sy n dv, Z Zz 
x x y y z z 
a, Oy 
Choose W% such that one of the — —— + 20v, = 0 and let the remaining two 
n, dv, df 
dv, be arbitrary, e.g., choose if such that pes + 2av_ = 0 and let dv_ and 
i ny, dv, x y 
4, Oy 
dv, be arbitrary differential variations. The implication is that a a 
Iv. 
| 


ov, = 0, i = x,y,z, and thus, nv) = Cyexp[-av,]. Finally, 


2 
-av : P 
= = iad dinates. 
n(v) nyt a fy ny) Ce in Cartesian vy-space coordinates 


Evaluation of the constant C -- 


fn) d2y here becomes 


mz 
u 


ma 
tt 


2 2 2 a, 3/2 
i i = (= T 
| { { Ce[-a(v, + vy + vy) dv, dv av, implies C CG) N 


=o co Sco 


Have used the gamma function defined by 


pe e" du, x > 0 
0 


T(x) 


and the properties 


T(x) = (x-1)P (x1) 


It 


T(m) = (m1)! m= integer > 1 


PGs) = vn 


The classical equilibrium distribution is therefore given by 


/2 L 


nv) = Go Ne 

Quasi-equilibrium classical distribution -- 

If n(x,v,t) is a “slowly varying" function of (r,t), then 
N(x,v,t) * ey N(r,t) era 


Recall that this result is in Cartesian v-space and thus might more appro- 


priately be written 


3/2 
= & ~o (v2 2 2 
A VeVyrV,.t) = C) N(z,t) exp[-a(v, + vy + v2)] 
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In spherical polar v-space coordinates, (v,0,), 
n'(v,0,6) = v? sing atv aviv) 


Whence, 


n'(r,V,0,¢5t) = @*? N(z,t)v2singe 


Particle speed distribution: 


n, (v) = [" ii n'(v,6) dédd 


0 0 


i.e., 


a 


ag? fl ¢2T 
@& i? Nv?e OV { { sin6dédo 
0 0 


a, (wv) 
Note that 


Tt f2T 2 
| [ sin@d6dp = [a Q = 4t 


1/2 3] 
a 


nj, (v) = an 


av 


Considering the relevant volume 
element for each case helps explain 


the different spectrum appearance. 


vos, ef 
Sor nd 
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Moments of the speed distribution -—- 


2 
ee Tae 


ven [ 
vrs vn, (vw)dv = 
y i (ra) 


|p 


E = ¢ f dgmv n, (v) dv 7 2 


Note that this relation will occur later in these notes as an identification 


of the constant a, i.e., 


ql 
eS [ee] 
bili 


Particle energy spectrum: 


' - dv 
n'(E) = n, [v(E)] dE 


3/2 1/2 = 
2 4 & Ne 20E/m 


Particle pressure: 


Rate at which particles in d3y 
strike a’s, as in the derivation 


of particle current relations, 


is given by n(v) dv viéds 


However, here assume that particles do not pass through d2s but are instead 


reflected. Presuming “specular reflection" (as illustrated), the particle 


momentum change in the surface collision is — 2mvcos0@, i.e., 


2mve& (-€). 


The magnitude of the force on d75 is 


2 
2m(vé) n(v)d *yd7s, 
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Thus, with p representing the pressure 


p= 2m(v*é) 2n(w) dv 
wre>d 


Case of an isotropic distribution -- 


p= 2m { (w°é)* n(v)d?v 
vre>0 


where n(v) is expressed in Cartesian v-coordinates. With the integration 


expressed in spherical, polar coordinates (@ as the polar axis), 


&S 
ut 


T/2 p27 
2m f { { (veos6)* n(v)v*sinddvdédd 


In particular, for the case of a classical equilibrium distribution 


2 3m, _ Nn 
pez eu GD * to 


Moreover, the "perfect gas law," or equation of state is p = NkT, where the 


Boltzmann constant is 


1.38 x 10746 


rs 
i 


erg/°K 


8.62 x 107? ev/°K 


and T is the “absolute temperature” of the particle distribution (in °K). 


Whence, 


ae 
2kT 


and, as a function of T, 


af2 wv? /2kT 
e m m 
nv, ¥¥5¥,) a NGue 


-mv* /2kT 


' ort a/2 2 
n'(v,8,6) “NG v’sin6e 


a 
a, (v) = WB Panton /2kT 
3/2 ae BA 1/2 
n' (8) = NGoE 4 WP hte Esker 


Elementary Quantum Statistics 


The classical equilibrium distribution is actually a limiting case (the 
classical limit) of more fundamental particle distributions. Because the 
energy states and energy exchanges in interactions for individual particles 
are of such small magnitude, the effects of the quantum (discrete) nature of 
energy exchange must be taken into account. In these lectures a heuristic 
derivation of such quantum effects and the resulting particle distributions 
are presented. 

Quantization rules: 

Wilson-Sommerfeld conditions —- 


Consider an independent set of coordinates x, which specify the configur- 


i 
ation of a system. If the forces acting on the system are conservative, then 
can define a set of generalized momenta p, by 
3 . 
Py =—77 Ex) . 


ox, 
i 


where E is the system kinetic energy and a = = X;° Moreover, if the sys— 
tem configuration is periodic (i.e., xy and Pp, are cyclical) then the "allowed" 
trajectories must satisfy a quantization rule fp za, = Oh where the integra~ 


tion is over a complete cycle of motion, a, is an integer, and Planck's 


i 
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constant is 


6.63 x yor? erg-sec 


a 
fl 


4.14 x 10" ev-sec. 


Examples -—- 


1. One-dimensional motion of particles between elastic walls. 


Implies p = mk 


9 
i oy 
Implies E on 


Periodic trajectory p = p(x) is 
as illustrated. The quantization 
implies that the shaded area (en-— 


closed by the trajectory) is ah, 


ines, 


2(2me)'/? x = ah 


Thus, the permitted kinetic energy states of the particle are given by 


242 
i week 


gmx 


a = 0,1,2,... 


Energy state separation is given by 


-~oa+i)ro  , dE 2oh? 
sux 8mX 


Note that the classical limit of continuously distributed allowed energy 


states is approached if either mor X is "large." Of more significance, 
= +* 0 as Q@ +, 


2. Three-dimensional motion of particles in an elastic box. 


x= (x,y.z) 
ke (x, ¥,2) 
E= jm (ik? +7 + 27) 


Implies P,, = mt, Ss = my 


Quantization rule implies that 


ah Qh ah 
Pe“ op? Py = “ope P2 > OE 


Whence, 
(of + ar 4 oy 
E pagename) Sabena 


A a Oo 
xyz 8mL 


Particle spin -—- 

In order for quantum mechanical arguments to yield meaningful predictions 
when applied to particle groups, it is required that certain particle types 
have an intrinsic angular momentum of fixed magnitude. This magnitude can 
have the values 


[s(s +1) GB*)/? 
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where s = 0, $ 1, 3, ... is termed the particle "spin." 
particle spin 
electron $ 
proton 5 
neutron $ 
photon 1 
(nuclei, ions, atoms molecules +, 3, 2, eee 


odd no. of elects., prots., neuts.) 


(nuclei, ions, atoms, molecules 
even no. of elects., prots., neuts.,) 0, 1, 2, 3, ... 

Spin states -- 

There is a "quantum number" associated with the spin orientation and 
this number can take on only 2s + 1 values (except for a photon where there 
are only 2 possible spin states with s = 1). An energy state is specified 
by trajectory quantum numbers and spin quantum numbers (e.g., for non-inter- 
acting particles in a box Os oy e are trajectory quantum numbers and for 


each set there are 2s + 1 possible spin states). 


Pauli exclusion principle ~- 
No two particles with half-odd~integer values of spin, s, can be in the 
same energy state (e.g., for non-interacting half~odd, integer spin particles 


in a box there are at most 2s + 1 particles with the same values of Oe oy? a). 


Thermodynamic probability: 
Particle distribution in energy -- 
Denote the permitted energy states of a system of particles by the set 


{E, i=1, 2, 3, ...} and the particle population corresponding to each 


i 
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state by {n, £02152 5: Byres nF 


—_|-Q— 
 ©- 

LQ 

—S-8-S- 


ji > 
£ = fe = - es 
: es &s fy “ 4, &@ % 
A= 7 A,tk& 4,°9 nye 7 4 = Yg= 
Heese Mg =O 
oh? 
Recalling the results for particles in a box of side L, Ey = Ser 


oh? -24 


and EY = mae For L = 1 cm., and m = 10 gram (ca. nucleon mass), 


AEY = 0(107 7) ev. Presuming that following energy exchanges of such 
small magnitude is not of interest in the particle distribution problem, a 
method of conglomerating individual energy states into groups of more signi~ 
ficant energy separation is considered. In the "macrostate"” description, 


the particle occupation population is the sum of the actual energy state 


occupation numbers, ny, over each conglomerate, e.g., 


—j—_~— 


| - 
inp # 
| i | nil 
/ v 20 28 
bine Bee Se E 
g, ey 
Ss 28 
Ht, = Pa Wy = 2s me 
a. &f=ZO 


The marcostate is described by the set {E,, N, i=1, 2, 3, ...}. Note 
that each macrostate is generated by many possible "microstates" which are 


described by the set {E,, n, i=l, 203 wesTS 


Ergodic hypothesis —~ 

Every microstate (consistent with any restraints on the particle system) 
is equally probable. Restraints on the particle system are such conditions 
as a fixed total number of particles, or a fixed total particle kinetic energy. 

Particle distribution in macrostates -- 

With every microstate presumed equally probable, the relative probability 
of finding the occupation number Ny in the macrostate ES is proportional to 
the number of microstate occupation possibilities consistent with (N,, E,). 

Let w(N,) represent this number of possibilities. The "thermodynamic prob- 


ability" of an occupation complexion is defined as 
N = Thy (N 
wt ep) Tw(N,). 


Role of restraints —- 
As an example, consider the case where the total particle population and 
total particle kinetic energy are fixed values. Then, certainly the range 


of possibilities for iN, Ej} is restricted by the constraints 


=z 
u 

mz 
ut 


om a 
i oa total population 


: Ny Ey = NE = total kinetic energy 
if there are further dynamic constraints, they must also be included. 
Expected particle distribution: 
The “expected" particle distribution is given by that occupation com-. 
plexion {N, } which maximized WIN, consistent with the relevant restraints 


on the system. Equivalent, and more convenient to apply, is the maximization 


of gn Wf ND . 
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In order to accomplish the determination of w(N,) and thereby wn, }), 
denote the number of microstate energy levels E;} in the ith nacrostate level 
(i.e., at E.) by M,. The counting procedure employed to evaluate w(N,) de- 
pends on the set of restraints on the particle system as well as a presumption 
on whether or not particles are distinguishable (i.e., in principle, can be 
labelled). Accordingly, different forms for the particle statistics emerge. 
In these lectures three forms are developed: 

1. Fermi-Dirac statistics. 

a. indistinguishable particles. 
b. half-odd-integer spin. 

ec. N fixed. 

d. NE fixed. 

2. Bose-Einstein statistics. 

a. indistinguishable particles. 
b. integer spin. 

c. N fixed. 

d. NE fixed. 

3. Boltzmann statistics. 

a, distinguishable particles. 
b. integer spin, or exclusion volume of particle occupancy. 
ce. N fixed. 


d. NE fixed. 


Fermi-Dirac statistics: 
A "Fermi gas" is a collection of particles with the properties 


1. they are indistinguishable (i.e., even in principle they 
cannot be labelled.) 


2. half-odd integer spin, 


3. fixed total population N, 
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4. fixed total kinetic energy NE, 


5. contained in a box of volume V. 


Note that E, here includés spin orientation, e.g., for s = 1/2, E, has 
E, and s "up," E, has E, and s "down." 
By the Pauli exclusion principle, the microstate occupation numbers, Des 


must be either 0 or 1. Thus, 
! 


and 
i! 
waN,}) =H Gap wt 
Whence, 


&n WEN, = % om, !-2n( M,N, )! - inN.t] 
and using Stirling's approximation (i.e., %nM! ~ M&nM-M for large M), 
== -(M, = = = 
anw({N}) CM, 2M (M N,)2n(M,. Ny) N20N 1 
which is to be maximized subject to the restraints 
IN, = N and IN.E, = NE 
yi £ ii 
Employing the method of Lagrange multipliers, maximize the function 


FC{N,}, y, &) = &n w({N,}) + y CRN -N) - B(EN, E,-NE) 


Note that 
aF M, 
aN, = 0 implies sae - 1) +y- BE, = 0 


if the approximation for mn WN, }) above is used. 


Whence, 


N 
i > lene (#Ey-)] +177 


Moreover, if this result is presumed to generate an approximate average occu- 


pation complexion for the My microstates, then 


n, * (exp{BE,-y] + 11> 


In quantum-number-space (i.e., G-space), 


ives 2s +1 
co exp [BE -YJ +1 


Where s is the particle spin and it is presumed that Ey is spin-orientation- 
independent. 
Let n(a)d7a represent the expected particle density in r-space and popu- 


lation in ao, volume element of a-space. Presuming that the particles are in 


a cubical box (of side L), 


n(a) = st 2s +1 
=“ 43 exp [BE -yltl 
Where 
. z z 
pe Qh _ ht 
2% 8m" ml.” 


Employing spherical polar coordinates in a-space (a,0,) 


m/2 ¢m/2 
nj (a) - { { n(a) a? sinededd 


a 0 
n, (a) = nar -y,, (Ba*h? -1 
al 13 (2s +1) f[e exp (E385 +1] 


Note that only the first octant of a-space is covered since 


a, = 0, 1, 2, ...; ay = 0, 1, 2, ...; a, = 0, 1, 2, ... 


In terms of particle kinetic energy and speed, 
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a" (E) = n, (a(B)) | $2] 

7 (2m3)1/2 gi /2 

= dn eo (2s + 1) tplke -y1 
ny (v) = al(E(v)) (= 


v 
ea” Yexp(B 4s mv") +1 


ft 


an®) *(s +1) 


Since n(Q) is "isotropic" in a-space, n(v) must be isotropic in v~space. 
oe O-sp: v, Ws 


Thus, in 6-dimensional Cartesian (r,v)~space, 


n(x,v,t) = ny (v) 


a 
4nv? 


@ Qs +1) (eYexp(8% mv’) +1)7+ 


It remains to identify the constants B and Y. The restraint conditions must 


be used, i.e., in the present terms, 


fn(z,v,t) dy = iE n, (v) dv =N 


and 
f{%s m’n(z,v,t) Ov =f 4auv n, (v) dv = NE 
Q 
The classical limit -- 
If for each macrostate N,<<M,, i.e., a very "dilute" occupation complexion, 
then 
! 
M,! 


vO) = GN DIN! 


M(M,~1) (M,-2). ar (M,-N, +1) 
N,! 


and M N 
ei i 
w¢tn, ) ; Wo 


L 


Presuming that N. is still large enough to use Stirling's approximation, 
N = IN M 
QnW({ 3) § 4 An! A faN, +1) 


and maximizing 2n W (14,1) subject to the restraints EN, = N and INE, = 
i i 


NE yields the result 


i = exp[Y-BE., ] 
L 
m, 3 7,78 le mv”, 
Whence, a, = exp[y-BE, 1, or n(r,v,t) = @ (2s +1) ele in Cartesian 
(x,v)-space. Note that this result is the classical equilibrium distribution, 


2 
i.e., n(x,v, » and that this classical limit is achieved by re- 


I< 
tt 
~~ 
g 
© 


quiring Ny << M,. Such a condition is developed if 


1. Nand L are fixed and E is increased, i.e., N, is decreased holding 


i 
M. fixed, or 
4. 
2. WN and E are fixed and Mor L is increased, i.e., M, is increased 


holding Ny fixed. 


Moreover, comparison with the classical equilibrium result yields the identi- 


fications 


and 


m3 7 m 3/2 
G) @st Ley = NGaer 


The identification of 8 is valid even for the case Ny = OM), however, the 


identification of y is only valid in the classical limit. 


Fermi energy: 


Defining W, the Fermi energy, by 


Mos YkT 


the Fermi-Dirac energy distribution is 


g/l? 
a'(—E) = Cc — 4 
EH) /kT 1 
3.2/2 
where ¢ = gr 222 (25 +1), 
h3 


KM = UCT) is determined from the restraint on particle populations, i.e., 


[ n!(E)dE = N 


9 


Consider the case of T > 0 and denote H(0) = Ho: The result is 


cK? for E< Ug 


n'(E) = 
0 for E>? Uy 

V2 

ate) CL : ae 
{ N= [ n'(E)dE = 5 Cuy 

| 

| Whence, 

| 

a 

oO Ay be cae gu? 


07 im “Gt@s Fi)! 


For the more general case, i.e., T # 0, numerical solution of N = [ n' (E)dE 


yields the results 
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Where u(T,) = 0 and kT) = 3.37 Uo 
For example -- 


Gas m N Yo T 
Molecular ct aes iol? 
(at STP) 


Electrons 1077? 1073 10 +5 
(in metal) 


i.e., for a molecular gas at STP, T > > Ty and - ~ >>] which imply that 
0 


a'(z) = col /kT,'4,-E/KT 


which is the classical equilibrium distribution. In contrast, for electrons 
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in a metal, T << T, and yp = uy which imply that n'(E) is like the T = 0 Fermi 


0 


distribution, 


For T = of 00°K) 
nile) f 


the electron distribution is 


EO as illustrated. 


Thermionic emission of electrons: 
Consider an electron gas in a metal with T = 0 (100°K). Suppose that 


the "wall" pee tree arsagnt the electrons in the "box" is such that an 


electron with. EB: > Ww hitting tie cecal muetace does escape from the metal 


volume. The "work function" of the metal is defined as d= W - Mo: 


z£ 


| i 
4 
_¥ 


tnside out s/he 


J, = ie ‘3 f- vip (vw) dv dv dv, 
2W. geo 00 , 
©) 


where 


n(v) = e" (28 +1) fexp[@mv*- y)/kT] + ryt 
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is the outgoing electron current at the metal surface. Using s =, U = us» 


: , 
and noting that ve implies Gmv? - 4) >> kT for ¢ = OCU), the result is 


2 
acy) = 2@)° expe ser] em /7Kt 


3 


2W. 
> (£8 
for vs ¢ =) 5 


whence, 


~ Atm 2 .~o/kT 
J + a (kT)“ e 


and the emitted electric current density is given by 


-$/kT 


I, = eJ, = 120 Te amp/em? 


4 
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where e = 1.5 x 10°” coulomb. This is the Richardson~Dushman relation for 


thermionic emission. 


Bose-Einstein statistics: 

A "Bose gas" is a collection of particles with the properties 

1. They are indistinguishable, 

2. integer spin, 

3. fixed population, N, 

4. fixed kinetic energy, NE, 

5. contained in box of volume V. 

There are no limits (other than the N and NE constraints) on the micro- 
state occupation numbers, nj. In order to determine w(N.) consider a sequence 


of Ni, periods (1) and NM ~ 1 commas (,), e.g., 


which represents the microstate complexion n, = 3, n, = 2, n, = 0, n, = 1, 


1 2 3 
ns = 3, . . . There are (N, + M, ~ 1)! permutations of the periods and commas 


of which Ny! permutations of the periods and CM, - 1)! permutations of the 


commas do not lead to different energy state occupations. Thus, 
(N, + M, - 1)! 
WN) = Qa Dp 
+ Nt CM, 1)! 


and 


(N, + M, - 1)! 
WIND IT a - DT 
i i 


Using Stirling's approximation, and maximizing {n WC{N.}) subject to the 


restraints 2 N. = N and 2 N.E, = NE yields 
. ~oRi 


i 
Ny 
a [exp(BE, - ¥) -1}7} 
Whence, 
a," (exp(BE, -y)- ay? 


and, with the identifications 


= 1/kT and y = u/kT, 
n, = fexp[(E, - Y)/kT] - iy? 


If the restraint 2 N, = N is removed (i.e., when particles can be created 
i 
or destroyed, then the parameter Y would not be introduced and the result 


is 


{expe /kT)- 2 


Photons in a cavity -- 

An electromagnetic wave of frequency Vv and velocity ¢ can be considered 
to be photon particles of energy hy and momentum hv/c with spin 1. There 
are two possible spin states associated with a auarod (mot the usual 2s + 1) 
corresponding to the two possible directions of wave polarization. The wave 


intensity is determined by the number of photons, and photons are easily created 
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and destroyed thereby generating a variable particle population problem. 
Consider photons in a unit volume cubical box with perfectly reflecting 


walls, i.e., 
Kare 2 2 2 2 
5) (a, + oy + on) h*/4 
or, 
2 2 o2(n2 2 2) 72 
(hy) c* (al + ay + a5) ne/4 
Ey = Rg = kohe 


Bose-Einstein statistics with a variable particle population apply, i.e., 


ny = 2 Lexp(hv,/kD) - ay 


Vy = “sc 
Whence, 
mo? 
a exp[hv, /kT] -1 
and 


da 
nj fa(v)] [SE 
vy 


n, () 


8mv2/c3 


ebv/kT 


The photon energy density is 


> s y 
From ob - ee” 
> 15¢e*h 


and the photon energy current per unit area is 


eo 5 4 
{ hyi, (v) cdv = 22D 
‘ 15c*h 
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A “black body" is an object which absorbs all incident radiant energy. Consider 
a black body placed in the box containing photons. When at equilibrium with 

the radiation in the box, the black body is at temperature T, and it is emit~ 
ting radiant energy at the rate of photon energy absorption, i.e., per unit 
area, the radiation rate is 


~ 2nek* 


ot” with st 
15c h 


which is the Stefan-Boltzmann radiation law and o is the Stefan constant. 


Boltzmann statistics: 

Consider a collection of N distinguishable particles in the sense that the 
labels 1, 2, 3,...,N can be meaningfully assigned to the particles (e.g., the 
atoms in a solid lattice). A microstate is no longer specified by merely the 
occupation numbers {n,}. Rather, it is now important to specify, in addition, 
where each labelled particle is in energy space. Presuming integer spin par- 
ticles, or exclusive volume of occupancy for each particle, there are My pos- 
sible microstate occupancies for each of the N, particles in the ith macrostate. 


Whence, 


Ny 
w(N,) = M, 


i 
Note that if WN, }) is taken as Fw (Ns), as used previously, then the inter- 
change of labelled particles between different macrostates is not counted as 
yielding separate and distinct particle occupation complexions. Consider a 
macrostate occupation complexion {N, } to be represented by a series of integers 


(the labelled particles) and commas (the separation between macrostate energy 


levels), e.g., 


1234,56,789,,4i10i11, 
For this case, N, = 4, with particles 1,2,3,4, Ny = 2, with particles 5,6, i = 3, 


with particles 7,8,9, N, =0... There are N! permutations of the integers, 
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However, permutations of the integers within a macrostate has already been in- 


cluded in w(N,). Thus, there are 


meaningful permutations of labelled particles with a given macrostate occupa~ 


tion complexion {n,t, and 


WON) = gay Twin) 
ede 
i.e, 
Me i 
wN,}) =N! q yi 


Maximizing eaw(tN,t) subject to the restraints EN, = N and IN,E, = NE yields 
i i 
Ne 
a, = exply-BE,] 


and 


a, = exply-BE, ] 


Identifying B = 1/kT and using aa = N gives the Boltzmann (quantum statistics) 
distribution 


exp[-E,/&T] 


i * N Sexp(-E, ke] 
L 


4-26 


Energy-Dependent Transport Relations 


As mentioned in the first paragraph of this section of these lecture 
notes, the E-dependence of particle distributions can be included in the 
formulation of the Boltzmann transport equation. So far in these notes, 
the derivations of relations describing particle transport were restricted 
to the presumption of monoenergetic radiation (pages 3-7 through 3-11). 
However, the arguments presented in the derivation of the monoenergetic 
equation (page 3-10) are easily extended to the E-dependent case with the 


changes 


v > v(E) 


(xr, 0,8, t) 


eo 
on 
nh 
re) 
. 
ct 
ww 
+ 


s(x &,t) * s(x,Q,E,t) 


ao 
¢ | areG'>+ ® + | | a7a'ae'e(8',E" > 4,5) 
4t 47 QO 


The emission distribution, f, has been redefined such that 


£(0',E" > 8,5) a70dE = the expected number of particles emitted 

in the increment a?Qan about (8, E) due to a collision involving 
a particle with pre-collision variables (',E"). Note that the 
total number of particles emitted per collision (c) in the mono- 
energetic version is now included in the requirements for f such 


that f is no longer normalized as a probability distribution. 
In these terms, the E-dependent Boltzmann transport equation is 


L a a “a a 
sty a (E,G,E,t) = - G-V0e,2E,t) - 2, (2.8, t)0(r,0,8,t) 


9 


+ | | ECE E',t)e(e,ts8",e" > Gz) 0,0" ,B,t)a7Mae" + 9(x,2,8,t) 


47 0 
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The solution complications introduced by describing the actual energy exchanges, 
which do, in fact, take place in most situations of interest, are certainly 
extensive. We shall approach the resolution of the E-dependent equation in a 
manner similar to that used for the monoenergetic relation, i.e., with simpli- 
fying presumptions which yield relatively easy to express notation and algebra. 
A conservation relation (the E-spectral equation): 
Consider the operation fan applied to the E-dependent transport equation. 


Term-by-term the results are 


aK 
v(E) dt 


O(x,B, t) 


12 FS 2 
TH) Be PEGE LER = 


| A+ 06 (c,0,£, t)d°Q = VeS(x,E, t) 


Z(E.E, t)O(,0,£, td = 2, (eB, tee, Et) 


foo) a Aa na 
| | ECE E' ,t)£(x,t3",n" > G,n)0(,8",B', tya"na7n" ae! 
4 


4n 47 0 


= [Petevas evetees8 > £E)O(r,E',t)dE" 
0 


| s(r,2,E,t)d-2 = s(x,E,t) 
where, 


O(r,E,t) = [owa.z.0a0 = the angular-integrated radiation 
flux at (r,E,t) which is also termed the flux E-spectrum 


at (x,t), 


JI(@,E,t) = [000.2 0)470 = the radiation current vector for 


particles of energy E at (r,t), 
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£(E' > E) = [ede > A,za’a = the collision transfer of energy 


description which has been presumed to be Q'~independent, 


s(x,E,t) = [samen = the angular-integrated radiation source 


density. 
Whence, the E-spectral equation 


a Wo a5 = 
VE) oe O CEE, t) Ve S(E,E,t) - 2, (2, E, to (x,t) 


+ |2,C2,E",t)f(x,t, E' + E)o(r,E',t)dE' + s(x,E,t) 
J 
0 


Thermodynamic equilibrium spectrum-- 

In the derivation of the classical equilibrium particle distribution 
(page 4-1), “equilibrium was defined by the restrictions: 

1) Steady state, i.e, t-independent 

2) Uniformity, i.e., r-independence 

3) Isotropic motion, i.e., Q-independence. 
To these we here add the presumptions of: 

4) No radiation absorption, i.e., ze = 0, or ae = ae 

5) No radiation sources, i.e., s = 0. 
The E~spectral equation becomes a relation describing thermodynamic equili- 


brium, i.e., 
EG)P) = [se yee > E)O(E')dE" 
(3) 
With reference to page 4-7 of these notes, we already have a solution 


of this equation, viz. 


@(B) = age =/*T 
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This should seem intriguing since we have not as yet specified 2 () and 
£(E' +E). The puzzle is understood by considering as a trial property for 


O(E) 


oe") EE! > E)O(E') = tte > E')O() 


This reciprocity (or symmetry) relation is termed the "principle of detailed 
balance." Its physical significance will be addressed in a moment. Here, 


note that substitution into the spectral equation yields 


£8) = 2,008) [ee + BtDa8" 
0 
Moreover, since the only allowed interaction is scattering, each collision 


results in the emission of a single particle, viz. 


fee > E')dE' = 1 
0) 
for this case. Thus, the principle of detailed balance is an acceptable 


property for the solution to possess. In fact, it is a required property. 
The physical meaning of the principle of detailed balance is easily 


illuminated by considering 
[2 "eB" ) de") [£@! > E)dE] 
= [2 (E)®(E) dE] [£E + E')dE'] 


which is obtained directly from the definition of the principle. The left- 
hand-side of the relation is a specification of the transfer rate of radi- 
ation from increment dE' at E' to dE at E. The right-hand~side is the 
transfer rate from dE at E to dE'at E'. The principle is thus a statement 


that in a thermodynamic equilibrium state each energy increment has a balance 
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of energy transfered to and received from each other energy increment. In 
this equilibrium state, the particles and targets must be distributed in 
energy such that the principle of detailed balance is maintained. This is, 
indeed, a property which the product 2, (E") £(E" +E) satisfies in the multi- 
scattering process which establishes the particle and target E-spectrums. 

It is important to realize that all the presumptions are required to 
yield a thermodynamic equilibrium condition. The particle moderation 
problem (considered next) is a good example of a situation where the equation 
to be resolved appears to be the equilibrium condition but results are clearly 
very different. 

Moderation spectrum-— 


A particle moderation problem is defined to be such that 
£(E' > E) = 0 for E > E' 


i.e., particles must lose energy in a collision. In order to achieve a steady 
state situation is a moderation problem, it is required to have radiation 
sources at some energy (usually high E) and radiation absorption at some 
energy (usually low E). To be specific, and to emphasize the comment made 

at the end of the discussion on equilibrium, consider the moderation case 


where all radiation sources have energy greater than E, and all radiation 


0 


absorption occurs at energies below EL with Ey < Eo: The E~spectral equation 


for a steady state, uniform, isotropic motion problem is identical to the 
thermodynamic equilibrium equation--with one important exception, i.e., 


E, <E< Ey is required. Specifically, 


2, (E)O(E) = [renee + E)O(E')dE" 
E 


for Ej <E< Ey 


where the lower limit of the integral reflects the moderation restriction 
applied to f. 

The moderation of fission neutrons by protons is a basic process in 
many thermal (spectrum) nuclear reactors. In this case Ey > 100 Kev and 
Ey < lev. The scattering interactions are elastic with isotropic emission 


in the collision center-of-mass coordinate system. These kinematic and 


emission distribution conditions yield 


£(E' +E) = 4, for 0<E<E! 


and the required £(E' > E) = 0 for E > E'. The moderation spectral equation 
becomes 
= v L ' t 
E(B) 2(B) exc ) Ez O(E')dE 
E 


<E< 
for E, E Eo 


which has the solution ®(E) = A/EX (CE) for E in the restricted range. This 
is not a classical equilibrium distribution and the principle of detailed 


balance is certainly not followed for E, < E < Ey (since only moderation 


1 
occurs). The equilibrium restriction of no sources or absorptions was 


relaxed and the solution is non-equilibrium. 


Plane symmetry problems: 

In order to pursue more non-equilibrium radiation transport problems 
and maintain a reasonably simple notation and algebra, consider restriction 
to plane symmetry problems in homogeneous, steady state media. The trans- 


port equation is 


1e 108 2 Gand 

VE) De POR HE,t) = — Ug OGeMLE,t) — 2 (B) OC, HE, t) 

° | f Z,(B')E(Q',E' > 9,8) O¢x,u" JE" ,t)d-R'dE" + 8 (x,U,E,t) 
4n 0 


wherein the integral term must clearly be reducible to a form which does 
not involve the direction components ¢ or ¢'. An often encountered case 


which leads directly to such a reduction is 
£0" zg! > Gz) = 2@-Q,e" > z) 


i.e., the case of azimuthally symmetric emission (about the pre-collision 
direction, 8"). Employing the usual form of Legendre polynomial expansion, 


i.e., 


2241 
4m 


£(Q'*Q,Et +E) = Y £,(E' + E)P, (A' +f) 
9=0 


and the spherical harmonics addition theorem, i.e., 
aie A <2 * 
1, ne ' ' 
PG") = apr 1 Ygg(t OX pq (ed) 


m=~2 


it is not difficult to show that 


CO 
| | ECE )EQ' AE! > E)6(x,u",E',t)d-Q"dE" 


4n 0 
co +1 

- | Be ew [ante ye, +2) | Pp (u") Gesu" VE", t) du! 
mo 0 i 


which demonstrates (for the case of emission azimuthal symmetry) the expected 
elimination of the @ and $¢' variables. 


The E~dependent, plane symmetry transport equation can thus be expressed 


as 
+ 9 o(x,u,8,t) = - pee O(x,u,E,t) - E,(E)0(x,u,£, t) 
V(E) Ot oH,B, Ug, Oke UE, t X,H,E,t 
S28 +1 
ye ey) | aBtz, (BD, (B" > B) { Po (ut )oGxsu' EB", t) du! 
0 ~1 


+ s(x,p,E,t) 
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and, the orthogonality of the Legendre polynomials can be employed to yield 


the E-spectral equation 


ae Mx,E,t) = - 32 T(E, t) ~ 2,(B) OB, t) 


v(E) ox 


+ f 2, CE") £q(E > E)O(x,E',t)dE' + s(x,E,t) 

It should be noted that the definition 
£(E' +E) = [ ect es >» AE)a22 
used in previous discussions of the E-spectral equation, and the Legendre poly— 
nomial expansion of the emission distribution yield the identification 
f(z >E) = £(E' > £) 

Whence, the principle of detailed balance takes the form 

D,CE')£,(E' +E) O(E') = L(E)£Q(E + E')O(E) 


These comments are included to indicate the particular physical significance 


of the emission expansion coefficient £)(e’ eB). 


Multigroup transport formulation: 
For the purpose of introducing the E-multigroup method of reducing the 
transport equation to more easily resolved terms, consider the E~spectral 


equation for a plane symmetry, steady state, uniform radiation problem, i.e., 


%, (E) &) = f 2, (E')£,(E! > E)O(E')dE' + s(E) 
i) 


Divide the energy range of interest into I contiguous intervals indexed 


as illustrated and define 
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ee 
v 
oP yo 
Eroup T aed Gr 6 
+ aes eee + Ss a Energy 
< £p7 Be. ey “o 
Ei 
o, = O(E)dE for i=1,2,..., I 
Ey 
The radiation with energy in the interval ES to Eo is termed group i radia- 
tion, and o. is termed the group i radiation flux. 
Bia 
The operation dE applied to the spectral equation yields 
E, E 
i ° 
Py = . 7 ' 
tifa Q(z > ij)O(E')dE' + S5 
0 
where, 
: es 
hae On | 2 (E)®(E) dE 
Ey 


is the group i total cross section (a flux-weighted average over the group 


energy interval), 


Bia 
Ey (Et + 4) = £,(B") fy! > B)dE 


Ey 


is the cross section for transfer (by radiation interaction caused emission) 


from energy E' to energy group i, and 


Fina 
8, = s(E)dE 


Ey 


is the radiation source emitted in group 1. To use consistent dependent vari~ 


ables, the integral term in the group spectral equation should be reexpressed 


according to 


a Lt 


E (B'+4)0(E')dE! = ») z (rie, 
0 


where the group-to-group transfer cross section, Er), is defined by 


Ey 


Ei) = ge | Dy (et+t) ec" ae" 
J 


Ej 


Note that these definitions are such that E,Gr4)9, is the transfer rate 
density of radiation with energy in group j to radiation with energy in 
group i due to interactions (not sources). 

If the group i radiation "removal" cross section is defined by ay = 
Dar - ZG), the usual form of the group spectral equation is obtained, 


i.e., 


Tri%i = J B Grade, + 8, 


for i =1,2,...,1 


This relation could be used to study E-spectrum problems and properties. 
Note, for example, that the group equivalent of the principle of detailed 


balance is 
I, G-) ‘> 2, (ind) o, 


if scattering is the only mode of radiation energy transfer and {4} is 
the thermodynamic equilibrium group spectrum. In these notes the deriva~- 
tion of the group spectral equation is intended to merely be an introduc- 
tion to the derivation and use of the more-complicated group transport 


equation. 


Multigroup transport equation-— 
Consider the E-dependent transport euqation for a plane symmetry, steady 


state problem, i.e., 


ge OGG ULE) +E, (EOC HED 


co) 


co 1 
2. 
= J 2 pay | al ce yeg(e're) | Py (u")OGx,ut, Bt) au" 
2% t L 2 
2=0 
0 -1 
+ s(x,H,E) 
Bia 
Applying the operation | dE, and the definitions of variables and 
Ey 


parameters used in the discussion of the multigroup spectral equation, yields 


the E-~multigroup transport equation, i.e., 


Wyk O Gu) +2, 8G) 


I +L 
5 Pou) Ey (ira) | Py (u')®, Cx," dh! + 8, (x,1) 
j=l aT 
where, direct extension of the definition of 2B Gra) is employed in the 


definition 


Biel Bia 
Ly (G74) = % de" 9(E")E, (E") £ ) (E'9E) dE 
E, E, 
The example of isotropic emission (i.e., z, Gea) = 2, Gri) 8,,) is somewhat 
academic since usually emission direction is correlated with energy exchange 
(through the kinematics relations). However, because of the resulting sim- 


plifications in algebraic manipulations, the isotropic emission case will be 


the example chosen (in these notes) for further amplification. 


If post-collision radiation emission is isotropic, the multigroup 


transport equation reduces to the form 


t tL 
a 1 , 
bar &G@.w +l, &G.w =F 2g (5>4) | 95% w)dul +s, (x,y) 


L < 
jel L 


for i = 1,2,...,1 


The structure of each member of this set of I equations is identical to that 
of the monoenergetic transport equation. The set is, however, coupled 
through the energy transfer term by the group transfer cross sections 
roid), and the monoenergetic flux, (x,y), is replaced by an "I-dimension 
vector" flux form {®, Ge, uw) 5 i=1,2,...,I}. 

In order to illustrate how the solution approaches, which were applied 
to the monoenergetic equation, can be directly extended to solution of the 
multigroup equation set, we shall here consider details of the Py approxima~ 
tion method. 

Solution by a Py approximation~~ 

Following the Pyvapproximation procedure as applied to the solution of 


the monoenergetic radiation transport equation, employ the expression 


N 
= 2nt+1 
2G, 4) = I Fr a PA 


in the I-group equation set, the relevant Legendre polynomial recurrence 


relation, and the linear independence of the Legendre polynomials to obtain 
nt 9 (x) + (ntl) Ho (xe) + (2M EL, 8 (x) 
dx “i,n-1 dx *i,n+1\* MEE) eg OE * 


I aor) 
= oh HGH C6, Hsp od 


for n = 0,1,2...,N and i = 1,2,...,1 


In the coupled set of (N+1)I linear, first-order differential equations 


S) n® is defined by a Legendre polynomial expansion of the source in 
3 


group i in a manner identical to the flux expansion coefficients, o, n® 
, 
and in equations indexed with n = N, o, nee is set equal to zero. 
> 


It is useful to reexpress the set of equations as 


d 


dx 1 Ba 4? ds 


9) - 1%, Ig (Grt)O, ox) = 5, 9 
I= 


and, for n > 0, 


pans) 


(x) + (2n41)2,, & |) =is. _ (x) 


(x) + (nt1) £ iow 


a 
ax o, n=l oe n+l 


The collision-induced radiation transfer-into~group-i term appears only in 
the n = 0 equations because of the presumption of isotropic emission. 

In order to continue this comparison of the E-multigroup solution 
approach and the monoenergetic technique, consider the homogeneous (s = 0) 
form of the transport equation with the purpose of determining eigenvalues 
and eigenfunctions. With the $5 ni) = 0 the equation set is translationally 


invariant with respect to x suggesting solutions of the form 


a -x/h 
a a) = ta 


The transport equation set becomes 


I 
¥ 1 - &, ¥ 0 + an LEGA, yO eS 


and, for n > 0, 


oY a) + CDE, 1 (© - (Qnty ez, HCO = 0 


for i = 1,2,...,I and n = 0,1,2,...,N with ¥ (2) = 0 in the n = N equation. 


i,N+1 


An example; two-group, diffusion analysis-— 


The example of I = 2 and N = 1 yields a two-group, P,-approximation 


1 
which is the simplest version of an energy-dependent diffusion approxima- 


tion. For this example there are four equations in the four unknowns 


Yo 62)s ¥y1(R)s Yog(2) and ¥y, (2). They are 


YD = lyr Yo (W) + Lg (272) ¥9q (2) = 0 


¥y 96%) - 3hE #440) =0 


i 
oO 


Yo, D - 8» Yo9 (2) + £2 (1-2) ¥, 9 () 


F596 - BRE 9 Yo) =0 


Note that the second and fourth equations constitute the two~-group version 
of Fick's Law of Diffusion. Non-trivial solution of the Yo n® implies the 
> 


secular equation which can be manipulated to the form 


4 2 
[z - Zy(L+2) 2, (2-1) 12 - {D,2.5 + Dt) +D,D, = 0 


cr? 12 


where d= 1/32) and Dy = 1/32, are the respective diffusion coefficients 


for group 1 and group 2 radiation. Whence, the four eigenvalues % = tL, and 


& = tL. where 


1/2 
— By 2 + Dore eeilts & 4D Do [24h 9729 (192) Ey (272) ] 
+ 212 Begg 172) By (271) J Dos + hen 


If either Eg (+2) =0or Ey (21) = 0, then it is more direct to rewrite the 


secular equation which can then be reduced to the form 


2 2 
(2k ~ Di) (Lk - D,) = 0 
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Whence, in these cases, the four eigenvalues 


= + = 
L t Ly and £= + Ly 


where 


are the respective "unclupled" diffusion lengths for group 1 and group 2 
radiation (and should not be confused with the same notation used for 
asymptotic eigenvalues in monoenergetic Py approximations, N = 1 and N = 2). 

Eigenfunctions and general solution forms vary significantly in expres- 
sion and physical interpretation according to the group transfer cross sec- 
tions. There are three cases and they are described here on an individual 
basis. 

Case 1: Ey (192) = (274) = 0. 

In this case the two energy groups are totally uncoupled and radiation 
transport is described by two monoenergetic relations. General flux solutions 


can be expressed as 


Yo 4 (4L,) 

U1 -x/L 

(x,y) = A [Ll + 34 pJe 1 
1 1 S10 CL ) 
Yo (-L,) 

ps cians & +x/L. 

+O +35 Ci.) HJe L 
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2 2 Vt) 
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+01 +3 yy ujet™/E 
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The eigenfunction ratios follow from the individual group expressions for 


Fick's Law, i.e., 


8 Cy 
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There is little need to further discuss the meaning of these solutions 


since they represent merely the superposition of two independent monoenergetic 


problems and thereby introduce no new physical concepts when compared to our 


previous monoenergetic radiation transport discussions. 


Case 2: 2) (12) # 0 and 2) (271) - 0. 


The opposite extreme relative to case 1, this case has complete coupling 


of the two radiation groups, i.e., transfer up in energy and transfer down 


in energy. 


® Cw) =A + (L+3qy 
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General flux solutions take on the more complicated form 


4-42 


As in case 1, the relevant eigenfunction ratios are determined by using the 


two-group form of Fick's Law, viz. 


Bice ae Ra © ha ae Oe 
hoGh) ore Hot Pe 
¥ 

Bg et ee Oy fait a 
¥ 5G) aes ice) = es 


Each energy group experiences influence from the other group (through 
the two nonzero transfer cross sections). Thus, the eigenvalues and eigen~- 
function ratios all carry parameter descriptions from both groups. There 
are no solutions which are strictly monoenergetic in character. 

It appears that there are too many arbitrary coefficients in the 
general solution form (Ayy> Cl A> Ci» Aga Coys Ay_> Cy_)- In analogy 
with monoenergetic diffusion results we expect two coefficients for each 
group allowing the satisfaction of two boundary conditions for each group. 
This yields a total of four coefficients rather than the apparent eight. 
However, substitution of the general solution form in the i = 1, n = 0 equation 
and in the i = 2, n = 0 equation gives relations between Ay and Ao.) Cay and 
Co. A. and Ay > and C,. and Cys Thus, there are actually available only 
the expected four arbitrary coefficients to meet boundary condition require~ 
ments for both groups. 

Case 3: U9 G2) # 0 but Ey (21) = 0, or Eg (2-1) # 0 but Zig (1+2) =0. 

The two possibilities lead to the same character of solution. Here, 
consider the usual relevant case of (2) # 0 but Eg (2+1) = 0. In this 
case, the i = 1 equations do not depend on the i = 2 solutions, but the 
i = 2 equations do depend on the i= 1 solutions. As in case 1, the eigen- 


values are Ly and Ly. General solutions take on the form 


4-43 
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The relevant eigenfunction ratios are again found by using the two-group Fick's 
Law. 

The group 1 radiation flux is independent of the group 2 flux (since 
%_ (2-1) = 0) and thus the solution form is the expected monoenergetic form 
with parameters descriptive of group 1 alone. The group 2 flux is dependent 
on the group 1 flux (since Eg (12) # 0) and thus the solution form is not mono- 
energetic in form and parameters descriptive of both groups appear. 

As in case 2, it appears that there are too many arbitrary coefficients 
(A. 


A, In order to satisfy boundary conditions we again 


gt Sie” Angi Pons Spy: Bag) 
expect only four coefficients. In this case, substitution of the general 


solution form in the i = 2, n = 0 equation yields relations between Ay and Mays 


and cy and Cor: Thus, again there are actually available only the expected 


four coefficients to meet boundary condition requirements for both groups. 
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ABSTRACT 


The mathematical formulation of a new approach to the homogéni- 
zation of certain types of heterogeneous media (such as a regular array 
of holes in a scattering media) for the purpose of neutron diffusion 
calculations is developed in this Report. The new method is based on 
the inclusion of an angular-dependent mean free path in the theory 
of neutron transport. In the present effort, calculations are restricted to 
media with plane symmetry and monoenergetic neutron theory is 
employed. 

It is found that a neutron-flux based theory and a collision-density 
based theory can lead to significantly different results when low-order 
approximations, such as diffusion theory, are employed in the solution 
of the transport equation. For the case of isotropic scattering, the 
normal mode technique is found applicable and closed-form, exact 
solutions are determined. 

There is little existing experimental work in the description of 
neutron distribution in such anisotropic media. It is therefore difficult 
to evaluate the results of the new theory. The interpretation of physical 
and numerical (Monte Carlo) experiments could be approached di- 
rectly with the new techniques. 


I. INTRODUCTION. 


Probably the most; important considerations of neutzon 


“This. Repoit..develdps: the: mathematical formulation -- 
diffusion in media s ‘holes are those é Behrens (Ref. 1). 


of a new sethad of tresting neutron diffusion in cértain’ 


types of heterogeneous ‘media. The heterogeneity ‘of in:- 


mediate concern, and toward which this work is directed, ~ 
is that cf a regular array of vacuum channels: (such as a .: 
* square lattice of cylindrical holes) in‘an otherwise homo- - - 
“geneous medium. The general procedure: with modifica?!" 

tion of- details, should be applicable to other. types. of-.: 


heterogeneity. However, 2 requirement which shouldbe 


- jimpdsed ic that the heterogeneity results in two charac: : 
‘s teristic directions, For example, in. the case of a regular ~ 


array .of vaccum channels, the two characteristic direc- 


* sNons- are parallel and perpendicular to the channel axis. 


“devoted to the determization of ths efeets.of heteragene--, 


» of the mediuia foz-neutron. diffusioet aloulations is. nat 


as the. process-wts: 


However,all previous work, including. thet in Ref 1; is 


ity on specific parameters relevant-te neuiron difusion, .. 
rather. than tc 2 formilstion of a general ‘method from 
which varisiss descriptiaris of the neutron, distyibutioiy can - - 


-be determined. Theigeneral problem nny be stated slog 3, 
these lines:is media, with beter: ageneity,. auch aa.the type 


considered here, it is clear that a “siraple’ bomegenization®, 


valid representation. A.simiple homeigwnizatien is, -defiped 
icinguthe medimadatoss sections’ 
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merely by the ratio of material volume to material-plus- 
vacuum volume. The streaming of neutrons in the vacuum 
channels leads to a spreading of the neutrons in the 
longitudinal direction (parallel to channel axis) which is 
larger than that in the transverse direction (perpendicular 
to channel axis). Thus, not only is the simple homogeniza- 
tion questionable for omnidirectional parameter calcu- 
lation, but also the anisotropic effects are completely 
subdued. Of course, the reason for considering homogeni- 
zation of the medium is the existence of an “arsenal” of 
possible mathematical attacks for such problems. 


The new approach to homogenization, which will be 
developed here, is based on the reasoning that neutrons 
traveling with a large component of their velocity in the 
longitudinal direction probably travel further between 
collisions, on the average, than those traveling with a 
large component of velocity in the-transverse direction. 
Let this effect be introduced into the neutron transport 
equation for the homogenized medium by allowing the 
mean free path to be angular-dependent. Clearly, by so 
doing, a mathematical fiction will be utilized since the 
mean free path, as used in neutron transport calculations, 
is a local parameter. This concept is certainly no more 
confusing than the idea of medium homogenization. It 
should be stressed that a total cross section which varies 
according to the direction of neutron travel, and not 


merely the usual scattering directional dependence, is 
being imposed. 


All considerations will be based on the idealization 
of monoenergetic neutron transport. Although time- 
dependent equations will be developed, most of the 
analysis will be devoted to the calculation of stationary 
states. For the purpose of completeness, the mathematical 
formulation will be developed along several lines. Spe- 
cifically: both the neutron flux and collison density will 
be considered as dependent variables; the familiar Py and 
double-P, approximations, as well as the moment decom- 
position will be applied; and it will be demonstrated that, 
for the case of isotropic scattering, the normal mode 
procedure, recently used for the solution of several types 
of neutron transport problems, is applicable and yields 
exact, closed-form solutions. 


The major part of this work is directed toward the 
mathematical formulation and physical interpretation of 
the new theory. A section, however, will be devoted to 
brief remarks relevant to application. There is meager 
experimental data available for lattices of the type con- 
sidered. It will therefore be difficult to evaluate the pres- 
ent theory. These considerations will yield a well-defined, 
albeit not well-substantiated, route to solution of prob- 
lems involving neutron diffusion in media pierced by 
vacuum channels. 


Il, MATHEMATICAL FORMULATION 


It would seem that application of these ideas to finite 

media dictates the use of non-separable position-angle- 

‘ dependent mean free paths. The inclusion of position- 

dependerice leads to gross difficulties which we have as 

yet not resolved. It will be assumed that the mean free 

- path in the homogenized medium depends only on the 

angle between the neutron velocity and the direction of 

the position variable. In all calculations plane symmetry 

is assumed, with the position variable either along the 
longitudinal or transverse direction. 


2 


A. The Neutren-Flux Equation 


The monoenergetic neutron transport equation for 
homogeneous media with plane symmetry may be written 
as ; ; 


1a 3 
O28 YOM) tas deat) + elu) y(en,2) 


=¢ foWF(2-2)y(xn',1) dQ’ +S(xut) (2) 


where y (x, x, t) is the neutron flux distribution as a func- 
tion of position x, direction cosine of neutron travel 
relative to x-direction », and time #; v is neutron speed; 
o(u) is the angular-dependent total cross section; c is the 
mean number of secondary neutrons which emanate from 
a neutron-nucleus collision; f (2° ’) is the normalized 
distribution in Q, the neutron post-collision direction of 
travel, of secondary neutrons produced by collision of 
a primary neutron with pre-collision direction 2’; and 
S (x, », t) is the rate of neutron introduction from sources 
which are independent of the neutron distribution. Al- 
though c and f are not necessarily descriptive of a non- 
multiplying medium, the terms scattering probability for 
¢ and scattering distribution for f will be used to avoid 
stilted discourse. An expansion of the scattering distribu- 
tion in terms of Legendre polynomials {P,, (Q+2’)} will 
be employed, i.e., 


faajy=ZElipiae) 


and the spherical harmonics addition theorem (Ref. 2, 
p. 148) to eliminate the azimuthal direction dependence 
appearing in the integral in Eq. (1). The result is 


19a C) 
Sam) + ase mt) + ou) ye mt) 


+1 
= FE Cat YhPals) [” PW ow oleae’ 
n bal 
+ S (x, u,t) (3) 
In order to proceed with a general discussion of the 
properties of Eq. (3) a more specific representation of 


o() is required. It is supposed, with little loss in relevant 
generality, that 


ols) = ZowPale (4) 

The sets {yn (x, t)} and (S, (x, ¢)} are further defined by 
trot)= ["Pae¥midn (a) 
Sa(t)= [Pale Slemt)dn (5b) 


Equations (5a and b) specify the respective expansion 
coefficients in Legendre polynomial expansions of the 
neutron Mux and source density, e.g. 


vlemt) = a+ on (2,1) Pa (n) (6) 
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In these terms, integration of Eq. (3) over the interval 
»é€(—1, +1) yields the relation 


Labo , ah &: 
0 at 7 ae tL) Z enn = So (7) 


Eq. (7) is the continuity equation for neutron motion. The 
only term which appears in an unfamiliar form is that 
which expresses the total neutron interaction rate. Clearly, 


-t 


E onda (2,2) = [ eerie (8) 


A further, familiar, reduction of the transport equation 
can be made in terms of the sets {yn}, {S,}, and (on). 
Using the recurrence relation for Legendre polynomials 
(Ref. 2, p. 32), 


(2n + 1) Pa(u) = mPa (a) + (A+ 1)Paar(n) (9) 


yields the set of coupled differential equations 
13 n 2a n+l a 
BaD TBE ag tO + ae Fe lon Ht) 
+ = (21 + )) omAimn ql = fn) yr (x, t) 
tm 


= S, (x, 2) (10) 
The set {Amn} is defined by 


Aime = 5 [Pile Palo) Pale de (aN) 


and has the following properties (Ref. 2, p. 87): 
1, The order of the indices is unimportant. 


2. Atma = 0 if the sum of any two of the indices is less 
than the third. , 


3. Aims = 0 if 1 + m+ nis odd. 

4, When Aim, 0, i.e., avoiding (2) and (3), 

Aes EY Otis) 
m= (To eaF) (2) (4) (+ m—n) 


(1) (3):--(+n~m-—1) 
x( Q@-Uta-m) ) 
(1) (3): (n+ m—1—-V) 
x( (Q) (4) (@tm—) ) 


(2) (4)---(L+ m+n) 
x (Geiser) (12) 
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It is from Eq. (10) that the various approximations to 
the transport equation are derived. Before proceeding 
with detailed examination of some approximations, re- 
examine Eq. (3). A simpler integral term is obtained if 
the dependent variable is changed to the neutron collision 
density 


F (x, ut) = o(u) y (x, 8). 


B. The Collision-Density Equation 


Reformulating Eq. (3) in terms of the collision density, 
F (x, », t), and the mean free path, A (u) = 1/o(), 


A(u) 2 
vo dt 


AF em) t+ awe F (x, n,2) 
+F (tpt) = 53 (Bn + 1) faPo (x) 


x qe Pa (u!) F (x, w’,t) du! + S (x44) 


(13) 

With the sets {F,, (x, ¢)} and {A,} defined by 
F(smt) = SA F* F(t) Pale) (a) 
A(s) = = AnPn (1) (14b) 


the set of coupled differential equations [cf., Eq. (10)] 
is obtained, 


= (21 + 1) Aim ae F, (x, #) 


+aaT “ T; = (21 + 1) Aim, wae Fi (x,t) 


n+1 
Qn+1 ; 


- (1 = fn) Fa (2, 1) = 


+ = (21 + 1) Aim,nsx Am =F (x, t) 


5, (x, t) (15) 


Note that if A, = 0 for n > 0, ie., the familiar case of an 


angle-independent mean free path, then Eq. (15), with the 


aid of the properties of {Aims}, reduces to 


Ao OF Mo OF a-1 , (n+ Dro OF as: 
o ot | ntl ot Qn+1 at 
+(L— ch) F, = 8, (16) 


which is the expected result. Note further that the n = 0 
member of Eq. (15) yields the continuity equation 


= So 
(17) 


Bt aa (21+ 1) wins EE + (1c) Fy 


In recognizing Eq. (17) as the continuity equation, use 
is made of the easily derived relations 


Yo (x,t) = 3S Ann (x,2) (18a) 


ts (% 4) = DY (22 +1) AmAnmEn (x, t) (18b) 


In the remainder of this work the case of a stationary 
state is assumed. In most cases this assumption merely 
leads to simpler algebra and notation and is actually not 
a requirement for the determination of a solution. The 
P,-approximation and double-Py-approximation, as ap- 
plied to both flux and collision density expansions, will be 
discussed and the moment decomposition for both flux 
and collision density will be considered. Finally, in greater 
depth, the case of isotropic scattering using the collision 
density as dependent variable will be covered. 


C. Fhe P;~Approximation 


The Py-approximation based on a flux expansion, or 
collision density expansion, is defined by the requirements 
that in Eq. (10), .or Eq. (15), ¥.(x) =0 for n>N, or 
F,, (x) = 0 for n > N, and the equations labelled by n > N 
are discarded. Thus, in a P,-approximation, N +1 
coupled differential equations are obtained with the 
N +1 dependent variables y, (x), n=0,1,-°-°, N, or 
F,, (x), n = 0,1, +--+ ,N. For example, the P,-approxima- 
tion based on the flux expansion gives the single relation 


oo(1 — c) Yo (x) = 0 (19) 


which has the implication c = 1 for non-trivial yo. This is 
a familiar result. The P,-approximation based on a col- 
lision density expansion gives the unusual relation 


aS i= (20) 


After this present brief comment considerations will be 
restricted to the case of A (x) a symmetric function of » on 
the interval » ¢€(—1, +1), and, in that case A, = 0. For the 
moment, suppose that A, 540. To be specific, suppose that 
4, >0. According to Eq. (20), the implication is that 


we DRM Ea ee 


dF ,/dx <0 which indicates a neutron flow in the +z- 
direction. The question which must be posed is: can F be 
independent (consistent with P,-approximation) and yet 
have ¢(x,) represent a neutron flow in +x-direction 
(ie. y increase with »)P Clearly the answer is in the 


affirmative since, if \(n) increases with » (implied by 


4, > 0), then the ratio y/A can be »-independent if y also 
increases with y. In all following considerations assume 
that A(u), and thus o(y), is a symmetric function of p. 
Thus it will always be required that 1, = 0, and ¢, = 0, 
for n odd. 


The P,-approximation (i.e., diffusion theory) based on a 
flux expansion gives the two equations 


Hs 5 (1 — c) ote (2) = 80) 


(21a) 
5G +0 -of)( 0 + Ze oe) = 5162) 
(21b) 


If it is further assumed that all neutron sources are 
isotropic such that S, (x) =0 for n>0, then Eqs. (2la 
and b) combine to give the usual diffusion theory 
relations 


dy 


— Doe + 0” ho (x) = So (x) 


on (22a) 


dbo 
dz —*(22b) 


ds (x) = 


with diffusion coefficient D, and “absorption” cross section 


o’, given by 
1 
D= 
3(1— of) ( o +-24) (23a) 
o = (l—c)ao (23b) 


It should be noted that only ¢ and o, enter into the 
diffusion theory parameters. The o (uz) expansion was not 
truncated at a quadratic in order to obtain Eqs. (23a and 
b). In passing, it is also worth noting, had time-depend- 
ence been included, that the P,-approximation would have 
given the “telegraphist’s equation.” The added assumption 
that dy, /vdt << 0y./dx results in the form of the familiar 
time-dependent diffusion theory with D and o’ again 
given by Eqs. (23a and b). 


- 1 
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The P,-approximation based on a collision density ex- 
pansion gives the two equations 


(24+ a) Ge + (0 — 0) Fale) = 8560) 
(24a) 

5 (d+ 3) GE +O of) = SG) 
(24b) 


In the case of isotropic sources, Eqs. (24a and b) combine 
to give 


2B 

(SE (2) = S02) (25a) 

= D’___ dF, 

Fi) = -— pv 

(2 + =) 
(25b) 
where the “diffusion coefficient” is 
(4238) 

D= (26) 


“3(L— chi) 


As is the case in the flux-based diffusion theory, a quad- 
ratic truncation of A(y) is not necessary to obtain the 
results given in Eqs. (25 and 26). 


The usual result of exponential spatial decline away 
from sources in infinite media is found for flux and col- 
lision density. The “period” of the exponential, the 
diffusion length L, is given by 


Li = Ty aN (27) 
Be (1 — 6) (1 ~ of) (ou + =B) 
based on a flux expansion, and 
(20 8) 
Li (28) 


*"30-o0—c) 


based on a collision density expansion. The dissimilar 
results obtained from a P,-approximation based on flux 
and collision density expansions have been pointed out. 
The divergence of results for the P,-approximation can be 
illustrated by considering the ratio Ly/Ly for a given 
problem. Suppose that A(4) is actually an even quad- 
ratic, ie., 4, = 0 for n = 1 and n > 2. The corresponding 
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(4yWte)® 


d2/ro 


Fig. 1. The ratio of diffusion lengths as calculated by a 
neutron flux based and collision density based 
P,-approximation for the case of an even 
quadratic mean free path, 

A Mw) = Ao + AgP2 {u} 


a(z) is not a quadratic, but only o, and o, need be com- 
puted since these coefficients determine Ly. In Fig. 1 the 
results are displayed for the range A;/A» € (0, 2). Clearly, 
the flux and collision density based expansions can lead 
to significantly different results. 


It should be expected that the two diffusion theories 
give different’ results. Note that whereas the y-approxi- 
mation yields an accurate representation for neutron 
current y, (x) but an inaccurate (truncated) representation 
for collision density, the F-approximation results in an 
accurate representation for total interaction rate Fy (x) 
but an inaccurate representation for current. For several 
reasons it will seem that the collision density expansion is 
favored in this work. However, these reasons are mainly 
of an algebraic character and it should not be construed 
that the F-formulation is, in all cases, superior. 


Approximations of higher order than P, are accom- 
plished following the usual general prescriptions. The 
added complications due to the angular-dependence of 
the mean free path place no restrictions on the formalism. 
Higher order approximations lessen the differences ex- 
hibited by the y and F-formulations. 


D. The Double-P,-Approximation 


- The double-P,-approximation is derived from Yvon’'s 
method whereby the angular-dependence of the neutron 
flux, or collision density, is decomposed into contributions 
from +2-directed and from ~x-directed neutrons. In 
order to simplify notation, let us consider the case of a 
stationary state in a medium characterized by isotropic 
scattering. The plane symmetry, monoenergetic neutron 
transport equation, Eq. (3), is then 


weve) +oW¥ ea =5 


x [ "a (a!) y (ea) du! + S (x, 2) 


(29) 
Use the half-angle-range expansions 
¥(on) = Z (2a + 1) ¥, (x) Pa (2p ~ Yn > 0 
(30a) 
= 5 (Qn + 1) ¥ (z) Pout Le <0 
(30b) 
S(x,n) = & (2n + I) Sy (x) Pa (Qu — 1), n> 0 
(80c) 
= E(2n +1) $; () Pp(Qe + 1),n <0 
(30d) 
o(u) = & 0% Pn (2x — 1), n> 0 (30e) 
= Do, Pa (2 +1),n<0 (30f) 
to obtain the set of coupled differential equations 
nave, ntldiin Wi 
Qn+ 1 de " Onti de “de 
+2 3S (21+ 1) Aton 9%, Hi 
ims 
= no S (oi vi + ov) + 25, 
(31) 


The double-P,-approximation is defined by the require- 
ment ¥%(x) = 0 for n > N and the equations labelled by 
n > N are discarded. 


The same analysis is followed for the collision density 
based approximation. Thus, with the definitions 


F (x, n) = & (Qn + 1) Fy (x) Pa(Qu— 1), >0 
; (32a) 
= & (2n + 1) F; (x) P, (Qn + 1),n <0 
(32b) 
A(u) = ZALPa Zn — 1),n > 0 (32c) 
= TARP a (Qe + 1), <0 (32a) 
the transport equation is obtained in the form 
2 aE; 
BET Zt Da Aimer GE 
n+l dFj 
+onF1 +1 E (2+ 1) As Amines Gee 
a SQl+VazA ail ai, op: 
= = im Aimn Fe 
= ¢8,, (Fi + F5) + 28% (33) 


The collision density double-Py-approximation is defined 
as in the y-formulation. 


In addition to the usefulness of Yvon’s method in the 
solution of problems where accurate representation of 
source or free boundaries are required, there is an added 
flexibility for the angular-dependence of the mean free 
path. Thus, A(z) may have one form for » > 0 deter- 
mined by the set {A*}, or {of}, and another form for 
» <0 following the set {Az}, or {o;}. For example, by 
using these methods a symmetric A(z) which varies 
linearly with » for » > 0 and » <0 can be expressed by 
using only two terms in each A(z) expansion, i.e., set 
Ags Ag, At = —Aj, and As, =O forn > 1. 


E, A Moment Decomposition 


In the usual theory of neutron transport through homo- 
geneous media, it is well-know that any space-angle 
moment of the neutron distribution can be found even 
though the distribution itself is unknown. In fact, an im- 
portant method of determining the neutron distribution 
is the construction of a likely flux shape from a finite set 
of moments. Let us now consider a moment decomposi- 
tion for the case of an angular-dependent mean free path 
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where, as in our previous considerations, two formulations 
will be examined, i.e., with y and F as dependent variable. 


We define the neutron flux and source moments by 


i= [elas (34a) 


5} = i ** IS, (2) dx (34b) 
Assume that the medium is of infinite extent, Multiplica- 
tion of the stationary state form of Eq. (10) by x/ and 
integration over x ¢(— ow, + ca) yields the set of algebraic 
moment relations 


B+ YowAinn (1 = ofa = Sh + BET sey 


x [ ays +(n+1) 2] (35) 


With the collision density moments similarly defined, 
ie., 


Fi= i "* IP, (x) de (36) 


and performing similar operations on the stationary state 
form of Eq. (15), the result is a set of algebraic equations 
relating the collision density moments, ie., 


A-cf) L=Si+ 5 


fod 
a Fi 


Pept hi} itn+1) S (+1) AwArm nea Ff? 


2n+1 fi, 

(37) 
The moments of the neutron distribution resulting from 
+a unit, plane, isotropic source (at x = 0) are easily in- 
terpretable in terms of important macroscopic parameters. 
For this source, S/,= 840840. Consider calculation of flux 
moments first and then contrast this result with the 
method applied to collision density moments. For the 
sake of definiteness, assume that ¢(z) is an even M de- 
gree polynomial in ». With o(n) an even function, it is 
clear that yi = 0 for 7 +n odd. An examination of Eq. 
(35), using the properties of the set {Aimn}, indicates that 
with M=40 and n+j even, y/ depends on 47%, $473, 
Waa Vie-wteas > Waya-es Youre Thus, yi,can not be 
determined without employing a truncation on the set 
{v»(x)}. This is, of course, the simplification used in a 
P,-approximation and would not yield exact moments. 
The familiar case of M = 0 poses no such difficulties and 

one can readily find the exact moments. 
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In considering the collision density moments, assume 
that A (u) is an even M degree polynomial in p. It follows 
that Fi = 0 for n + j odd. It is also possible to show, from 
Eq. (37) and the properties of {Aim}, that F/=0 for 
n>j(M +1) and therefore, with finite M, any collision 
density moment can be calculated exactly. For example, 
for the case of M = 2 


E 18 4g 
ce ) a 
<p ttre -Jja-qy | 
where (x*) is the mean value of x? and Lr is the diffusion 
length as calculated by a collision density based 
P,-approximation. Note that the result of Eq. (38) is unlike 
that found in the angle-independent mean free path case. 
In that case M = 0 and the second spatial moment is given 
correctly by diffusion theory, i.e., (x?) = 2L* both in the 
exact calculation and in diffusion theory. 


In passing, note that the set {yi} can be found from 
the set {F/)} via the easily derived relation 


w= = (21 + 1) AmAtmn FY, (39) 
The result of Eq. (39) does not contradict the earlier 
assertion regarding the problem of finding the flux 
moments. When the determination of the {y/} is ap- 
proached directly by Eq. (35), it is the total cross section 
which is considered to be an M-degree polynomial in p. 
When {y/} is found using {Fi} as in Eq. (39), the mean 
free path is assumed to be an M-degree polynomial in p. 


F. The Case of Isotropic Scattering 

The stationary state form of Eq. (13) with the additional 
assumption of isotropic scattering gives the collision 
density equation 


(a) oF (en) + Fu) = $f. Fewhdw +50) 
(40) 


In this case of isotropic scattering an angle variable 
change is suggested. Specifically, define the angle vari- 
able u = pA(u)/A(1), measure x in units of (1), and 
change dependent variable to F (x,u) = g(u) F (x, p) 
with g(u) =|dp/du}. With the source density change 
S (x, u) = g(u) S (x, »), Eq. (40) takes the form 


ut F(s,u) + F(x, u) = ew [ Peuydu + S(x,u) 


(41) 


Various aspects of the solution of Eq. (41) will be 
considered. 


1. Legendre Polynomial Expansion 


Following the procedure used in the derivation of Eq. 
(10), the coupled differential equation form of the trans- 
port equation 


n 4@Fy-y +i CLavee 
Sit a esr a tee 
_ bn 
=e F, (x) + S, (x) (42) 


is obtained where the Legendre polynomial expansions 
are used: 


F (zu) = SE* rey Pa(u) (4a) 
S(u)= SAS* s(2)P.(u) (43) 
g(u) = B-esPa(u) (48e) 


It should be noted that the requirement of a symmetric 
4 (u) imposes the condition that g, = 0 for n odd. 


The idea of a Py-approximation is equally well-applied 
here. For example, the P,-approximation (diffusion theory) 
takes the usual form [cf., Eq. (25)] 


1 dF, 
a 3 a + ad- C80) Fy = So (44a) 
_ _ 1 dF, 
oS ae 


where only isotropic sources are allowed. 


2. Moment Decomposition 

With the collision density and source moments defined 
in the usual manner [i.e., by Eq. (34)], Eq. (42) can be 
transformed to the algebraic set 


CEn 
3 = S} ————— Fi 
Fi = Si+ Sai 
(45) 
+ Brey nFis + n+ DPE 


It has been assumed that A () is symmetric which implies 
that g(u) is symmetric. We also find that Fi=0 for 


j +n odd, and, for odd n, F{ depends only on Fi-* and 
Fi;1. Furthermore there is the interesting property that 
the spatial moment FJ depends only on the set of moments 
{F!, n + isS7}. Therefore, the calculations of a low-order 
spatial moment require the specification of a small num- 


ber of the g, and the prior determination of a small. 


number of other moments. 


As an example, calculate by these methods the second 
spatial moment of the neutron distribution resulting from 
a unit, plane, isotropic source (at x = 0). In this case 
$4 = 80540. The moments F? F®, and F? are easily deter- 
mined and are the only values required in the calcula- 
tion of F2. For the normalized second spatial moment 
Icf., Eq. (38)], 


2g, 
Gets (46) 
3 3 1L—cgo 


3. Normal Mode Expansion 


The recently developed normal mode technique (Ref. 
3) will be applied to the problem of determining the 
exact and asymptotic solution of Eq. (41). In so doing 
there is generated an interesting mathematical problem 
the details of which are considered in Appendix B. Con- 
sider the homogeneous form of Eq. (41), ie, S=0. 
Translational invariance suggests the “ansatz” 


F (2,1) = 6 (1) exp (—) (47) 


where the separation variable v is allowed to be complex. 
The integral equation 


(=u) 86,4) =F vg(u) a $(v,u")du’ (48) 


‘is thereby obtained. Adopting the usual normalization, 


f rte ee (49) 


If solutions of Eq. (48) are allowed to be distributions (in 
the sense of L. Schwartz, Ref. 4), then 


$44) = £78 s aaw—u) (60) 


Assume that g(u) satisfies a Holder condition (Ref. 5, 
p. 11) on the interval of the real line ue(—1, +1). Any 
singular integrals which might appear are then of the 
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Cauchy type and their evaluation is defined as the Cauchy 
principal value (Ref. 5, p. 26). 


The normalization required of ¢(v,u), ie, Eq. (49), 
leads to a specification of allowed discrete values of v in 
the region of the v-complex-plane excluding the line 
(~—1, +1), and to a specification of the function A (v) for 
ve(—1, +1). To aid in the analysis of these results, define 
the Cauchy integral, G(v), by 


GW)=55 i) é fie. rT (51) 


With v¢(—1, +1), Eq, (49) gives 
1 + irevG (v) = 0 (52) 


which has a set of roots which are distinct. With 
ve{—1, +1), Eq. (49) yields an explicit formula for the 
function A (v) and no restrictions are placed on allowed 
values of v. We find 


A(v) = 1+ inevG (v) (53) 


Thus, if the definition of A(v) is extended, as expressed 
in Eq. (53), to the entire v-plane, the zeroes of A (v) deter- 
mine the set of allowed distinct values of v. Since g (u) 
is symmetric, G(—v) = —G(v), whence, A(v) is an even 
function of v. The zeroes of A(v), therefore, appear in 
pairs which are labeled =v;. 


A set of functions of the angle variable u indexed by v, 
{$(v, u)} have been found. There is a discrete indexed set 


with v¢(—1, +1) and members characterized by 


— © ¥18 (4) 

$ (by; 4) = Qyse (54) 
and, a continuous indexed set with ve(—1, +1) and of 
form given by Eq. (50). The function A (v), which appears 
in ¢(v,u) for ve(—I, +1), is given by Eq. (53). Further- 
more, the zeroes of A (v) for v¢(—1, +1) establish the set 
of discrete indices {-:v;)}. 


If we assume that g(u)=40 for ue(—1, +1), Eq. (48), 
with the normalization of Eq. (49), may be written in the 


form 
[:-= = (58) 
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Multiply Eq. (55) with index v by ¢(»’,u) and sub- 
tract the result from Eq. (55) with index v’ multiplied 
by ¢(v,u). Employing Eq. (49) and integrating over 
ue(—I, +1) 


[-s1) ae o(v,u)e(v’,u)du=0 (56) 


There is clearly no degeneracy and thus Eq. (56) may be 
rewritten as the orthogonality relation 


[eee du=Oforvsey” (67 


The nature of the orthogonality relation including the 
case v =v’ depends on whether v is a member of the 
discrete index set or belongs to the continuum. If v is 
a discrete index, then 


= Faw PlRYHt) a(evy w)du = 1 (sy) 855 (68a) 
2, [7 _ug(u) 
sa daar is (vj u)* Gee <AGE) 


If v belongs to the continuum, then 


“ou A = vA? (v) lth 
i PC) Rah ll aaa ZG) 8(v— Vv’) (59) 


It was found that the set of normal modes {¢ (v, )} is 
orthogonal, with weight function u/g (u), on the interval 
ue(~1, +1). For the remainder of this section assume 
that the norma] modes are also complete in the space of 
functions which satisfy a Holder condition on the interval 
te(~1, +1). Appendix B, in measure, substantiates this 
hypothesis by demonstrating the existence of the modal 
expansion coefficients. In so doing, the interval of com- 
pleteness is generalized to include all physically relevant 
cases, 


Assuming that the normal modes form a complete set 
on the interval ue(—1, +1), the general solution of Eq. 
(41) is in the form 


F(qu)= Sav)e(nuer(—=) (60) 


where the summation indicates integration over continu- 
ous spectra when applicable. In many problems there are 
boundary conditions which can be formulated as 


F (0, u) 


= @(u) = ZS a(v) }(v,u) forue(—1, +1) 


(61) 
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and the orthogonality relations can be used to determine 
the expansion coefficients, a(v). In detail, Eq. (61) is re- 
written as 

¢(u) = Ealtys) (+054) + > a(—vj) 6 (—4,4) 


+ [7 a0) 6 (ow) deforwe(—1, +1) (62) 


Direct application of the discrete index orthogonality 
relation, Eq. (58), yields the discrete indexed expansion 
coefficients, 


1 a ou 
a(ztv5) = Teal ata)  (u) > (sev, u) du 


(63) 
Using Eq. (57) in Eq. (62) gives the result 
+1 7a = +1 u 
[FeHweeomae= [an ooo 
x i a (v’) cy (v’, u) dy’ 
(64) 


There appears a doubly singular Cauchy integral and thus 
the order of integration in Eq. (64) may not be reversed 
without due caution. The doubly singular term appears as 


c +1 

vad 3 ; du——— 
Assume that a@(v) satisfies a Holder condition for ve(~1, 
+1) and follow the dictates of the Poincaré-Bertrand 


formula for inverting the integration order (Ref. 5, p. 57). 
Then 


ug(u) [7 vale 
y~ Use, vm 


dv’ 


+1 +h: al , 
of dy 28) wel) ay = oF vg (v)a(v) 


au ZO) [HB a, 65) 


v—u “1 saa 


2 
+5 i‘ 
Using Eqs. (59) and (65) results in the more useful 
“orthogonality relation,” 


+2 u +1 ; ‘ ; 
[ dur 604) [ a) eblwar 


= I(v)a(v) forve(—1, +1) (86a) 
romeo [(38) *(8)] om 


Now, applying Eq. (66) to the problem of finding the 
continuum expansion coefficients in Eq. (62), 


1 +1 tu 
a) =f sey eeb.edde — 6D 


for ye(—1, +1) 


G. The Green's Function for the Case 
of Isotropic Scattering 


As a specific example of the use of the relations just 
developed, consider the problem of finding the infinite 
medium Green's function for isotropic, plane sources. In 
this case, the source density, S (x, u), of Eq. (41) represents 
a unit, plane, isotropic emission of neutrons at a position 
which is designated x=0, ie, S$ (x,u) = g(u) 8 (x)/2. 
Integration of Eq. (41) over a vanishing interval about 
x = Q yields the boundary condition 


© 


u[F (0*,u) — F (0-,u)] = forue(—1, +1) 


(68) 
The additional condition that as [|x| 0, F(x,u)—> 0 is 
imposed and the solution is expressed in the form 


F (s.u) = E.a(+n) 6(+ynu) exp(—*) 
i id 


ee a (+) $(+,u) exp (—* *) drfors > 0 
(69a) 


F (z,u) = — 2a(—v)$(—vp u) exp(=) 


= a (0) $(,u) exp (—) dv forx <0 
(69b) 


The source condition, Eq. (68), then takes the form of the 
general boundary condition, Eq. (62). Specifically 


Zep Tate (tou) + 7 Tw al— 9 (— yn) 


glu, 


u a _itl 
taal a) $(,u)de= 2 
(70) 
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Whence, employing the normalization expressed by 
Eg. (49), 
1 


a(=v)) = We) (718) 
atv) = ary forve(- 1,+) (71b) 


and the solution of the Green’s function has been com- 
pleted. 


Some aspects of the Green’s function will be examined. 
For simplicity, assume that there is only one pair of zeroes 
of A(v), stv. A sufficient condition for this property is 
developed in Appendix A. The Green's function is then 
given by 


Fee) Soe e (=) 
; aa. exp (=) dvforx > 0 
° (72a) 
F (au) = ~ $7228) exp (=) 
: ee exp (=) dy forz <0 
i (72b) 
With the definition 
én(v) = 2 $ (1,4) Py (u) du (73) 


and the easily derived symmetry properties, 


I (vo) = —I(+vo) and I(—v) = ~I(v), 


the collision density moments for the neutron distribution 
from a unit, plane, isotropic source, are 


ofl ” 
Fi = -y [Gn (tv0) + (— 1) on (— Nay 5) 


+ FL te) + (De FG ae 
4) 
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Using Eqs. (9), (43c), (48), and (49), a recurrence relation 
is obtained for the set {$x (v)}; 
(2n + 1) udu (v) = ndn-a (v) + (+ 1) guar (v) + cv, 
(75) 


and the normalization ¢»(v) = 1. Note that Eq. (75) im- 
plies that ¢,(v) is an even, or odd, polynomial in v of 
degree n. Therefore, ¢,(—v) = (~1)"¢n(+v), and Eq. 
(74) reduces to 


FL = all py en(to) + { “Fey eae | 


ifj+niseven (76a) 


ifj+nisodd  (76b) 

Already considered is the moments set {F4}. For this 
particular case the {F/} is determined from Eq. (45) and 
the source condition Si = 8y8j. In passing, note that the 
consistency of Eq. (45) and (78) is easily demonstrated via 


the recurrence relation Eq. (75). Moreover, equating the 


F° and F3 moments as derived by the two relations gives 


Tw ak Tor! "tae (tia) 
i i+ “oes 
T(¥ 0) +[F TO) 74 =3 (1 — cgo)? (77b) 


From Eqs. (77a and b) an explicit expression for the dis- 
crete index v, is obtained: 


NGhie a 


_ 3 G—ca)? (l- ra 


= ae Toy 4 


For c<1, vo is real and is interpreted as the exact 
asymptotic diffusion length [here, measured in units of 
4(1)]. It should be noted that the integral terms in Eq. 
(78) depend on ¢ and {g,} via the dependence of I (v) on 
these parameters [cf., Eq. (66b)]. 


(78) 


Il, REMARKS REGARDING APPLICATION OF THE THEORY 


A limited number of considerations relevant to the 
application of the theory presented in Section II will be 
developed as a brief illustration of possible methods of 
application of the present theory to physical problems. 
Many interesting calculations are possible, and with the 
accomplishment of experimental measurements of neutron 
distributions in the types of media under discussion, many 
comparisons of theoretical and experimental results would 
be profitable. 


Methods of determining the proper variation of mean 
free path are certainly required if these mathematical 
formulations are to be applied to physical problems. In 
this section the general types of heterogeneity, toward 
which the current theory applies, will be discussed. A 
simple method will be detailed, using known diffusion 
lengths, to specify the angular dependence of the mean 
free path for a particular type of heterogeneity. 
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A. Types of Heterogeneity 


The motivation of the present effort, as mentioned 
earlier, is the establishment of a method of homogeniza- 
tion of regular arrays of vacuum channels for the purpose 
of neutron diffusion calculations. Also imposed was the 
necessary restriction that, in general, the type of hetero- 
geneity considered should yield two characteristic 
orthogonal directions. As an example of the caution which 
must be exercised in application of the theory, consider 
a type of heterogeneity which, at first approach, appears 
to satisfy the necessary requirements, but which actually 
is unsuitable for these methods. Specifically, examine the 
ease of a periodic slab array of scatterer and vacuum. 
This heterogeneity exhibits two characteristic orthogonal 
directions; perpendicular to slab, and the directions in the 
plane of the slab. Moreover, the direction perpendicular 
to the slabs (transverse to slab “channels”) yields con- 
siderations which are algebraically easily accomplished. 


If A(u,x) represents the mean distance traveled to a 
collision by a neutron located at a position x to the left 
of the right-hand-face of a slab of scatterer, traveling with 
direction cosine , relative to the slab perpendicular direc- 
tion, gives 


forxST,,2>0 
(79) 


In Eq. (79), A, is the mean free path in the scatterer 
material which has slab thickness T,, and T, is the 
vacuum slab thickness. For the homogenized medium we 
require a function A(z) which, it would seem, should be 
a “suitable” average of (n,x). For the case of isotropic 
scattering, the average 


"2 (an) Yo () de 


A(u) = ae 
f yeloide (60) 


is clearly indicated. In Eq. (80), 1 (x) represents the actual 
angular integrated neutron flux. Note that, in the present 
case, Yo (x) can be found. Far from neutron sources there 
is Yo (x) ~> exp (x/L,) where L, is the “asymptotic” diffusion 
length in the scatterer material. A note in passing: if 
T,/liag << 1, then yo (x) is approximately constant and 
Eq. (80) gives the result 


T, 


aw) =a(1+ =) (81) 


which is the “simply homogenized” parameter. Of course, 
the condition T,/L, << 1 should yield the homogeneous 
limit. 


If x now represents the direction perpendicular to the 
slabs we have the asymptotic result y. (x) > exp (~x/L,) 
when the position x falls in a scatterer slab and y (x) 
is a constant when x falls in a vacuum slab. The “best 
fit” to this flux, for the homogenized medium is 
o (x) > exp (—x/L) where L is the simply homogenized 
diffusion length, ie, L=L,(1+T,/T,). This result 
would be obtained if Eq. (81) were used. In this partic- 
ular case, the situation is that a calculation based on an 
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angular-dependent mean free path yields results less rep- 
resentative than the simply homogenized calculation. It 
is expected that, in the orthogonal characteristic direction 
(i.e., in the plane of the slab), use of an angular-dependent 
mean free path is indicated. 


The case of a calculation in the slab perpendicular 
direction for a periodic slab array is certainly excluded 
from the present considerations. Moreover, no motivation 
should be felt toward developing a theory for that case 
since it is easily treated by a standard method, i.e., change 
of position variable to “optical thickness.” Consider now 
the details of a macroscopic-parameter-based calculation 
for a heterogeneity for which the present methods were 
clearly intended, i-e., a regular array of cylindrical vacuum 
channels. 


B. Cylindrical Channels in a Regular Array 


Consider a regular array of vacuum channels of cylin- 
drical cross section. With every vacuum channel of cross 
sectional area A, a cross sectional area of scatterer mate- 
rial A, is associated such that V = A,/A, is the ratio of 
vacuum volume to scatterer volume characteristic of the 
medium. Label the axial, or longitudinal, direction with 
x and direction cosine », and the radial, or transverse 
direction with y and direction cosine 7. Due to streaming 
along channels, different diffusion properties in the x- and 
y-directions are expected with both of these cases being 
different than the simply homogenized diffusion. The 
simply homogenized mean free path is given by 


ds = A,(1L + V) (82) 


Before presenting a specific method for obtaining a repre- 
sentative A(u), some general comments will be made 
regarding the features of such a calculation. It is clear that 
only two representative directions are being considered: 
axial, or x-direction; and transverse, or y direction. It is 
also clear that in the present lattice the actual description 
of a straight line path in the transverse direction. starting 
from a point in the scattering material depends upon 
both the azimuthal angle about the x-direction and the 
particular position in the scattering material relative to, 
say, the center of a vacuum channel. A “suitable” aver- 
aging technique must be employed. Furthermore, the 
same difficulty is encountered when considering a descrip- 
tion of an axially-oriented path, Let A. (4) and A, (9) rep- 
resent the angular-dependent mean free paths with 
respect to x-direction diffusion and y-direction diffusion. 
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The following constraints on the “suitable” averaging 
technique seem intuitively reasonable: 


1. The average mean free path based on axial and 
transverse directions should be equal, ie., 


f * heli) da = f Ba (a) dy (83) 


2. The axial mean free path in the transverse direction 
(i.e., at ~ = 0) should be equal to the transverse mean 
free path in the transverse direction (i.e., at 7 = +1), 


ie, 
Az (0) = Ay (£1) (84) 
3. Both axial and transverse mean free paths should be 
symmetric, i.e., . 
Ae (nu) = Ae(—») (85a) 
Ay (9) = dn (—7) (85b) 


A possible method of obtaining A(z) is to find, as a 
function of starting position in the scatterer, the mean 
free path length traveled in all directions. Then, upon 
“suitably” weighting this quantity, according to whether 
Az(u) or Ay{n) is desired, an average yields the angular- 
dependent mean free path. In even the simplest lattice 
this is a geometric task of considerable magnitude. Here, 
for the sake of brevity, an alternate, albeit certainly less 
self-contained, route will be taken, Assume certain macro- 
scopic diffusion parameters, such as diffusion length, are 
given and use the general constraints of Eqs. (83-85) to 
obtain a representation of the mean free path which yields 
the given parameters. To be specific, assume that 4, (z) 
and Ay (y) are even quadratics of the respective variables. 
Thus, in terms of the Legendre polynomial expansion 


Ye (1) = Azo + AzzPo (u) (86a) 


Ay (m) = Apo + AysP (x) (86b) 
From Eq. (83) Azo = Ayo, and this result used in Eq. (84) 
yields Ayz = —Azz/2. Therefore, in terms of the two un- 
knowns, A. and A,, Eq. (86) may be reformulated as 


Ae (nu) = Ao + A2P2 (yu) (87a) 


y(n) = de =-5 MPs (0) (87) 
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The arguments used here with respect to the mean free 
path also apply to the determination of the total cross 
section. Thus, the general constraints are expected: 


1 f . oe(n) dp = £ “ (n) dy 
2. oe, (0) = oy (2:1) 


3. oe (p) = oe(—p) 


and 
ay () = oy(—y) 


If the total cross section is assumed to be an even 
quadratic, then in terms of the two unknowns, o» and o:, 


Ge (n) = ao + oP 2 (n) 


1 
y(n) = 0 ~ “5 oP s(n) 


For the remainder of this discussion assume that the 
neutron collision density is used as dependent variable 
and thus the mean free path is the relevant parameter. 
These considerations can be equally well applied to the 
neutron flux and total cross section. 


From Eqs. (28) and (87) results are 


e 5 


ie Is (88a) 
~~ 3(L— ¢)(1 — cf) 
(.-#) 
Peer iieaael i ae (88b) 
" 3(L—¢)(1 — of) 
Also, 
peeks eee (89) 
* 3(1 =e) (1 = ef) 
From Egs. (88 and 89) 
do _ (L) 
qa. (80a) 
(L) = = + “ 7 (90b) 
de 57 L, Ly 
=-3[2-F| (90e) 


Measured values of L, and L,, or other theoretical treat- 
ments, can be used to find L, and L, in order to determine 


Ao/As and As/A, For example, if Behren’s theoretical 
formulation (Ref. 1) is used, 


(#) =1+2v+—2 __ + anv 
exp (F =i : 
(91a) 
(=) = P+ av + ca +RV 
ow (7) -1 
(91b) 


' where R is the ratio of the vacuum channel radius to A,. 
In Fig. 2 L/L, and L,/L, are presented as a function of 
R, Re(0,5), for the cases V = 0.5, 1.0 and 2.0 as deter- 
mined by Eq. (91). Then in Fig. 3 for the same values of 


Fig. 2. Neutron diffusion lengths as found 
by Behren’s theoretical 
formulation 
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R and V, the results for 4./A, and A,/A, are based on the 
curves in Fig. 2. 


It should be noted that what is referred to as L? and L3 
is Eq. (91) are actually calculated by Behrens (Ref. 1) as 
(x*)/2 and {y*)/2 and, via Eq. (38), 


2) = of: +18. 
=i a= jana) = 
AR 
8 4 
Wy) = 2p +e (92) 


175 (l—c)(l— cfs) 


If AZ > 1 the validity of Fig. 3 as a relevant representation 
for (u) is questionable. However, truncation of 4 () at 
a quadratic would, in that case, also be of questionable 
usefulness. 


Fig. 3. Expansion coefficients for an even quadratic 
representation of the mean free path 
angular dependence 
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IV. SUMMARY 


The mathematical formulation has been developed of a 
new approach to the homogenization of certain types of 
heterogeneous media (such as a regular array of vacuum 
channels) for the purpose of neutron diffusion calculations. 


The new method is based on the inclusion of an angular- 
dependent mean free path in the theory of neutron trans- 
port. In the present effort, calculations are restricted to 
media with plane symmetry and monoenergetic neutron 
theory is employed. Extension to energy-dependent theory 
and to other symmetries would probably follow the gen- 
eral lines for the familiar, angular-independent case 
without significant additional complication. However, it 
seems clear that the requirement of the existence of two 
orthogonal characteristic directions in the development 


of the angular dependence of the mean free path must 
be imposed. 


It was found in this Report that a neutron flux based 
theory and a collision density based theory can lead to 
significantly different results when low-order approxima- 
tions, such as diffusion theory, are employed in the solu- 
tion of the transport equation. For the case of isotropic 
scattering, the normal mode technique is applicable, and 
exact, closed-form solutions can be determined. 


Evaluating the results implied by the present theory 
with respect to measurements is impossible. There is a 
current lack of pertinent experimental results for neutron 
distribution description in the relevant type of media. 


APPENDIX A 


The Function A (v) 


It has been found previously that the zeroes of A (v) for 
v¢(~1, +1) are the discrete set of normal mode indices 
and that they appear in pairs, +v,;. The number of these 
allowed discrete indices will be discussed. To this end, 
and for relationships which are useful in Section II G, 
there follows a brief study of the general properties of the 
function A (v) as defined by Eqs. (51) and (53), ie., 


A(v) = 1 + inevG (v) 


“ “6 a 
me fi, wy 


Gi) = 


In terms of the set {g,}, as defined in Eq. (43c), A (v) may 
be rewritten 


A(v) = 1— cv = EnQn (v) (A-1) 


where ©, (v) is a Legendre function of the second kind 
defined for the entire v-plane by an extension of the 
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Neumann formula (Ref. 2, p. 51) to include v¢(~1, +1), 
ie., 


Qn(v) = + if 7 En) du (A-2) 


1 ¥ 


with singular integrals evaluated as the Cauchy principal 
value. For large v, vQn(v) varies as v-", Thus, A(v) is 
bounded for large v. Furthermore, the Q, (v) are analytic 
in the v-plane excluding v e(—1, +1) and, therefore, A (v) 
is analytic in this same region. The contour illustrated in 
Fig. A-1 and the argument theorem (Ref. 6, p. 116) estab- 
lishes the number of zeroes of A({v) in the region 
v¢(—l, +1). Since the zeroes of A (v) appear in pairs, the 
number of zeroes are denoted by 2J. The argument 
theorem applied here yields 


4nJ = change in arg A*(u) on C, 


+ change in arg A- (u) on C_ (A-3) 


-i Ct. | 
| 


Fig. A-1. Contour in v-plane used in determination of the 
number of zeroes of the function A (v) 


It was assumed that g (tu) satisfies a Holder condition on 
uwe(—1, +1) and therefore G(v) is a Cauchy integral. 
The Plemelj formulae (Ref. 5, p. 43) are applied to find the 
limit values G* (u). It is found that 


G(u) = G(u) = > alu) (A-4) 


where G*(u) and G(u) refer to the limit values of G (v) 
as v approaches u from above and below the real line 
respectively. From Eq. (A-4) the limit values 


At (u) = A(u) = “cug (u) (A5) 


Now, A(u) with we(—1, +1) is a real function (with 
singularities at u = +1), so that A(0) = 1 and g(u) is a 
symmetric function. Whence, the relations 


arg A*(u) = — arg A~(u) (A-6a) 
arg A*(0) = arg A-(0) =0 (A-6b) 
At (u) = Av (—u) (A-6c) 


These results used in Eq. (A-3) yield the number of pairs 
of zeroes J in terms of the single angle arg A*(+1), ie., 


J= 4 arg A*(+1) (A-7) 


It should be noted that Eq, (A-3) contains the implicit 
requirement that A*(u)3<0 for ue(—1, +1). This as- 
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sumption is not completely necessary; however, it prob- 
ably applies to most cases of physical interest and its 
application greatly simplifies these considerations. 


A sufficient condition for J = 1 will be developed in the 


~ case that g(u) is an N-degree polynomial in u, i-e., 


glu) = & exPa(u) (A-8) 


Note that the Legendre functions, Q,(v), can be ex- 
pressed as 


Qn (v) = Pa(v) Qo(v) — Wa-r(v) (A-9a) 
Qo(v) = arctanh vy, ve(—1, +1) 
= arctanh—>,v4(—1, +1) — (A-9b) 


where W,.-1(v) is an even, or odd, polynomial in » of 
degree n — 1 (Ref. 2, p. 51). In these terms A(v) is re- 
written as 


A(v) =1—cvQo(v) & gaPn(v) + a gnWa-s (v) 
nad ns (A-10) 


Also, as u> +1, Qo(u)—> +00 and P, (+1) = 1. Clearly, 
W,,.(+1) is bounded, and thus, if 


3 g.>0 (A-11) 


nao 


then as u> +1, A(v)-> — o. From Eq. (A-5), in the 
present case, 


At(u) = A(u) + Fen 5 GnPn (tt) (A-12) 


Therefore, the following may be concluded: if Eq. (A-11) 
holds and, in the range u e(0, +1), 


3S gaPs(u) >0 (A-13) 


then arg A*(+1) =a and we have the desired result, 
J = 1. It should be stressed that Eqs. (A-11) and (A-13) 
give a sufficient, not necessary, condition for the number 
of pairs of discrete indexed normal modes to be unity. 
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APPENDIX B 


A Relevant Hilbert Problem 


In Section II F the existence of the modal expansion 
coefficients {a (+v;), j= 1,2, °-+ ,J,a(v), ve(—1, +1)} 
was assumed. Moreover, the orthogonality relations are 
based on the whole angle range ue(—1, +1) and thus 
only provide a means of determining expansion coeffici- 
ents for the case of a boundary condition given over all 
angles. It is found, by reducing the problem of finding 


expansion coefficients to the solution of an inhomogeneous 
Hilbert problem, that the existence of expansion coeffi- 
cients can be demonstrated, and a method prescribed for 
determining the value of the coefficients for problems in- 
volving all physically relevant boundary conditions. The 
techniques elegantly described by Muskhelishvili (Ref, 5), 
are followed closely. 


I. REDUCTION OF TRANSPORT PROBLEM TO AN INHOMOGENEOUS HILBERT PROBLEM 


Transport problem boundary conditions will be en- 
countered, in general, of the form 


$(u) = = a(+nja(Pna+ z a(—vy,1) $(—»5,4) 


+f a(v) } (v, u) dv for u € (a, 8) (B-1) 


where —la< @<=+1. If it were possible by some 
method to determine the set of discrete indexed coeffi- 
cients {a (2:v;),7 = 1,2, ++ - ,J}, and define 


# (uw) = 60) — 3 altn)a(toyu) 
— 3 a(—w) 6(-vou) (B-2) 


then there would be an integral equation for a(v), 
ve (a, 8)» Le, 


[P eWermar=swforuetae) 83) 
Using the derived form of ¢ (v, u) Eq. (50) 
A(u)a(u) +$ew [Oe 
= ¢’ (u) for ue (a, 8) (B-4) 


From Eq. (A-5) 


Suga) = sy lar(w)— aw] Ba) 


A(u) = [aru + a-(u)] (BS) 


and therefore Eq. (B-4) may be rewritten as 


[a (u) + A°(u)] ua (u) + [ar(u) ~ 4° (A) 
= ud’ (u) (B-6a) 


for u€(a, 8) 


There is assurance that A(u), as defined in Eq. (B-6b), 
exists if a (u) satisfies a Holder condition on u €(«, 8). For 
the moment, assume that this condition is fulfilled and 
define the Cauchy integral, A (v), over the entire v-plane, 


(B-6b) 


8 
al f? welt) oy 
Qari J, uy 


Aly) = (B-7) 
The Plemelj formulae yield the limit relations on the line 
ue (a, 8 ), 


At (u) ~ Av (u) = ua (u) (B-8a) 


At (u) + Ar (u) = 2A(u) (B-8b) 


The results of Eqs. (B-5) and (B-8) applied to Eq. (B-6) 
give the alternate form 
At (u) At (u) 


— A (u) A (u) = u¢’ (u) for ue (a, 8) 
(B-9) 


It was assumed that A*(u)=<0 for we(a, 8). With this 
restriction Eq. (B-9) can be easily transformed to the 
form of a boundary condition for an inhomogeneous 
Hilbert problem on an arc (Ref. 5, chapter 10). The prob- 
lem of determining a(u) restated in these terms is: 
Find the sectionally analytic function, A (v), vanishing at 
infinity, with boundary condition on the line ue(a, 8), 


A-(u) 
a (uy 


ue = 
(u) 


At (u) = = A-(u) + (B-10) 


Note that the assumptions on g(u) and A*(u) imply that 
A (u)/A*(u) is a function satisfying a Holder condition 
and not vanishing on ue (a, 8), and, if it is assumed that 
the angle boundary condition, ¢(u), satisfies a Holder 
condition and a (=v;) exist, then ud’ (u)/A* (u) satisfies a 
Holder condition on u €(«, 8). 


Summarizing will help clarify the procedure. If it is 
assumed (what is to be proved) that a(u) satisfies a 
Holder condition on ue(a,f), then the integral A(v), 
defined by Eq. (B-7), is of the Cauchy type. Now, Cauchy 


JPL TECHNICAL REPORT NO. 32-686 


integrals are sectionally analytic functions with boundary 
the line of integration. Specifically, if (a, 8) is the line of 
integration: 

1. A(v) is analytic in v-plane excluding (e, 8). 


2. A(v) approaches well-defined limits as we (a, 8) is 
approached from either side of (a, 8) with possible 
exception of the end points, « and/or 8. 


3. Near the end points, A (v) satisfies the conditions 


JAW Soe as 
[AW <r asy 


where a, b, A and B are real constants, and a <1 
and b <1. 


Moreover, A(v) vanishes as |v|-> oo. The integral equa- 
tion for a(u) has been transformed into the boundary con- 
dition Eq. (B-10) which is the form of an inhomogeneous 
Hilbert problem boundary condition. Thus, the original 
transport problem has been reduced to an inhomogeneous 
Hilbert problem. If a solution A (v) which introduces no 
physical ambiguity can be found, then the assumption of 
the existence of a(u), ué(a, 8), will be substantiated. 


Hl. SOLUTION OF THE HILBERT PROBLEM 


In terms of 
8 (u) = arg At(u), A> (u)/A* (u) = exp [—2i0 (u)] 


the Hilbert problem boundary condition [cf., Eq. (B-10)] 

becomes 

of (u) 
At (u) 


At (u) = exp [—26 (u)] Ar (u) + for u € (a, 8) 


(B-11) 


Since A (v) must also vanish|v|-> co, the solution (Ref. 5) 
is 


Hv) eee) du 


AW) = “a7 "= Gye wre) 


(B-12) 


where H (v) is the fundamental solution of the associated 


homogeneous Hilbert problem and is given by 


H (v) = (a — vj 8/* (B — y)@B/* eo) (B-13) 
The Cauchy integral @ (v) is defined by 
vane Eee) 
Ov) = =f rer du (B-14) 


Providing x = 6 (8)/m — @ (a)/z is a positive integer, there 
are the x additional requirements 


8B ubrig’ (u) = = re _ 
} (a) Hu) du=0 forn=0,1, wold 


(B-15) 
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These additional requirements are a necessary fea- 
ture of the solution. It should be recalled that the 
function ¢’(u), ue(a,f), is not completely specified, 
i.e, the discrete indexed expansion coefficients, 
a(+v;), in Eq. (B-2) are, as yet, unknown. For 


Il, APPLICATION OF THE 


Plane symmetry transport problems fall into two gen- 
* eral categories: 


1. Infinite media problems with full-angle-range bound- 
ary conditions (such as the Green's function solved 
in Section II G). 


2. Half-space media problems with half-angle-range 
boundary conditions (such as albedo or Milne type 
problems). 


Combinations of the solutions of these type problems 
lead to the solution of cases with finite media (slabs). For 
full-range boundary conditions, the orthogonality of the 
normal modes provides a direct method for determining 
expansion coefficients. The solution of the Hilbert prob- 
lem in these cases demonstrates the existence of the co- 
efficients and thus partially supports the completeness 
hypothesis. For half-range boundary problems, there are 
no apparent orthogonality properties of the normal modes. 
In these cases, the solution of the Hilbert problem not 
only provides proof of existence, but also gives a well- 
defined prescription for the determination of expansion 
coefficients. The application of the Hilbert problem solu- 
tion will now be outlined to the categories of full-range 
and half-range boundary conditions. 


In the case of an infinite medium, full-range boundary 
condition problem, a source condition is usually given at 
some position, which is chosen to be designated x = 0. For 
e <1, it follows that F (x,u) should vanish as |x|-> 0. 
Thus, the general form of solution is that given in Eq. (69). 
The source condition can be formulated as 


pu) = 3 altn) o(tenu) + 3 o(-n)o(—enu) 


+ [arenas forue(—1, +1) 
: (B-16) 
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the general problems considered later, it will be 
demonstrated that, in each case, the x require- 
ments are necessary and sufficient for the complete 
specification of all discrete and continuum expansion 
coefficients. 


HILBERT PROBLEM SOLUTION 


Instead of using the obviously indicated orthogonality 
properties, consider the coefficient evaluation by the route 
prescribed in the Hilbert problem solution. Note that 
a=—l and B= +1. From Eqs. (A-6) and (A-7), the 
results are @(—1) = —Jm and 6 (+1) = Jx. Therefore, in 
this case, x = 2] and there are 2/ requirements of the form 
of Eq. (B-15). Specifically, 


+1 untig’ (u) os 
a aa) ay 2 = 9 


Equation (B-17) provides a sufficient number of equations 
to find the discrete indexed expansion coefficients, 
a(+v;), 7 =1,2,°-° ,J. The fundamental solution H(v), 
{cf., Eq. (B-13)], is given by 


forn=0,1,--°°,2J-1 
(B-17) 


H() =(-1L-vy 1—»e®™ = (B-18a) 


O(v) = me ean du 


ol, aay (B-18b) 
Thus, A(v) is determined [by Eq. (B-12)] and @(u) for 
ue(—1, +1) can be found from the limit relation [ef., 
Eq. (B-8a)] 

ua(u) = At(u)-— Av(u) — forue(—1, +1) (B-19) 
Since the problem has been completely and unambig- 
uously solved, it is clear that the supposition that @(u) 
satisfies a Holder condition is substantiated and the 
existence of the expansion coefficients has been dem- 
onstrated. 


For half-space media, consider two types of problems. 
An “albedo problem” is described by a boundary con- 


dition at the medium surface (x = 0 with medium occupy- 
ing x > 0) specified for u¢ (0, +1) and the condition that 
F (x,u) vanish as x «. A “Milne problem” is described 
by a similar boundary condition at += 0, but with 
F (x, u) > $ (—v, u) exp (x/v) with v = v,,7=1,2,---,J, 


or ve(0,+1), as x-» 00. These problems have been 


specified as boundary conditions on the half-range 
ue(0, +1), With obvious modifications, the procedure is 
easily applied to half-space media occupying x < 0 and 
boundary conditions on ué(—1,0). With the half-space 
occupying x >0, the general solution of an albedo 
problem is 


F(x,u) = 2 a(+v) 6 (+, u) ex (+) 


+ i a(v) p(v, u) exp(—*) dv forx > 0 
(B-20) 


and for a Milne problem, 


F(x,u) = A$(-»u)exp(+) 


v5 


+ ‘ "alo #(nu)exp (—*) do forz > 0 
(B-21) 
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In both cases, the boundary condition at x = 0 can be 
expressed in the form of Eq. (B-1), ie., 


blu) = E alte oronu) + fab) slaurar 


for ue (0, + 1) (B-22) 
Now, a = 0 and 8 = +1 and, from Eqs. (A-6) and (A-7), 
8 (0) = 0 and 6(+1) = Jz. Thus, x =J and we have the 
J requirements 


+1 urrig’ (u) = = . < 
f eres hc aa forn = 0,1, Jol 
(B-23) 


These are sufficient to determine the discrete indexed 
coefficients, a(+v;), j=1,2,--- ,J. The fundamental 
solution takes the form 


H (v) = (1 — v) exp (@ (v)) (B-24a) 
eee -=f . Lo) ay (B-24b) 


The Hilbert problem solution, A (v), v¢(0, +1), and the 
continuum expansion coefficients, a(u), we(0, +1), are 
found as in the case of a full-range boundary problem. 
Again, substantiation is found for the supposition of the 
existence of the relevant members of {a(v)}. Moreover, 
a prescription is found for calculating the expansion 
coefficients when the use of orthogonality conditions is 
impossible. 


NOMENCLATURE 


@ modal expansion coefficients of F 
A Cauchy integral with density a 


A, cross-sectional area of scatterer material per lat- 
tice cell 


A, cross-sectional area of vacuum per Iattice cell 
Atnn 
¢ neutron scattering probability 

D diffusion coefficient 


triple Legendre polynomial integral 


f neutron scattering distribution 

fo Legendre polynomial expansion coefficients of f 

F neutron collision density 
F, Legendre polynomial expansion coefficients of F 
Fi angle-space-moments of F 


g Jacobian of angle-variable change, i.e. 
&(u) = |du/du| 
g. Legendre polynomial expansion coefficients of. g 


ai 
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NOMENCLATURE (Cont'd) 


Cauchy integral with density g 
fundamental solution of Hilbert problem 
orthogonality integral value 


half of the number of discrete indexed eigen- 
modes 


neutron diffusion length based on an F-expansion 
neutron diffusion length based on a y-expansion 
asymptotic diffusion length in scatterer material 
diffusion length with respect to x-direction 
diffusion length with respect to y-direction 
Legendre polynomial ; 

Legendre functions of the second kind 

ratio of vacuum channel radius to A, 

neutron source distribution 

Legendre polynomial expansion coefficients of S 
angle-space-moments of S 

time variable 

slab thickness of scatterer material 

slab thickness of vacuum material 

angle variable defined by u = pA (u)/d (1) 
neutron speed 

ratio of vacuum volume to scatterer volume, i.e. 
V = A,/As 

regular part of Legendre functions of the second 
kind 

position variable 

mean value of x? 

position variable perpendicular to x-direction 
Dirac 8-distribution 


spt rr Fr >» 


vn 
vn 


Kronecker 5-distribution 

included in the interval 

direction cosine of neutron travel relative to 
y-direction 

argument of the function A 

Cauchy integral with density 6 

index of the Hilbert problem 

neutron mean free path 

simply homogenized medium mean free path 
Legendre polynomial expansion coefficients of 
mean free path with respect to x-direction 

mean free path with respect to y-direction 
coefficient of $-distribution in angle eigenmode 
direction cosine of neutron travel relative to 
x-direction 

angular flux eigenvalues 

discrete indexed angular flux eigenvalues 

total cross section 

absorption cross section 

Legendre polynomial expansion coefficients of o 
total cross section with respect to x-direction 
total cross section with respect to y-direction 


angular flux eigenmodes (and angular boundary 
conditions) 


Legendre polynomial expansion coefficients of 
neutron flux distribution 

Legendre polynomial expansion coefficients of y 
angle-space-moments of ¥ 


unit vector in direction of neutron travel 
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